OpenGNSSLabSPP
Flagship tutorialNumerical walkthroughRINEX + Broadcast EphemerisReal data examples

SPP from Scratch

Compute the GNSS receiver position from raw pseudorange observations, step by step.

Level

Intermediate

Estimated time

2-3 hours

Parts

Part 1 of 3

Selected SV

G09

GNSS satellites orbiting Earth illustration

Geometry snapshot

G09

Satellites in view

5

PDOP

3.15

What is Single Point Positioning (SPP)?

Single Point Positioning estimates a receiver's coordinates and clock bias from code pseudoranges and satellite orbit/clock information. This worked GPS example follows real observations through transmission time, broadcast propagation, corrections, linearization and iterative weighted least squares.

Prepared by OpenGNSSLabUpdated Sources, reproducibility and corrections

Tutorial overview

Numbers before abstraction

This tutorial builds a Single Point Positioning (SPP) solution from first principles using only code pseudorange observations and broadcast ephemeris. Each formula is paired with real numbers so you can see exactly how results are obtained.

  • Derive the measurement model and apply key corrections
  • Linearize and solve the system with weighted least squares
  • Iterate to obtain the receiver position and clock bias
  • Validate the solution and understand common pitfalls

What you build

A clear, modular SPP pipeline you can implement in Python, MATLAB, or C++.

Processing pipeline

From code ranges to receiver position

RINEX observations

Raw pseudorange measurements

Satellite position & clock

Broadcast ephemeris + clock model

Earth rotation correction

Sagnac effect applied

Atmospheric corrections

Ionosphere & troposphere

Linearized model

Observation equation

Weighted least squares

Solve for state update

Receiver position

Iterate to convergence

Selected inputs

Raw observations and selected satellite

RINEX snippet

Raw observation epoch

Selected satellite: G09
SVPseudorange (m)ElevationAzimuthSNR
G0920,149,603.51274.43 deg55.08 degnot exported
G0721,134,422.50465.60 deg307.44 degnot exported
G0421,325,674.30147.10 deg99.20 degnot exported
G3022,518,702.19936.90 deg268.55 degnot exported
G1623,464,421.18019.60 deg42.58 degnot exported

Broadcast ephemeris

G09 compact row

Toe

378,000 s

t_tx

377,999.932345 s

sat clock

4.430247e-4 s

more parameters

NAV record

Elevation

74.43 deg

Azimuth

55.08 deg

Strong signal / good geometry

Checkpoint 1

Worked numerical arc (G09) - Pseudorange model

Pi=ρi+c(δtr−δti)+Ti+Ii+di+ϵiP_i = \rho_i + c(\delta t_r-\delta t_i)+T_i+I_i+d_i+\epsilon_i
Measured pseudorange20,149,603.512 m
Geometric range20,282,408.337 m
Satellite clock term132,815.465 m
Sagnac correction-5.640 m
Ionosphere9.999 m
Troposphere0.000 m
Speed of light299,792,458 m/s
P_corr = 20,149,603.512 - c*dt_s - Sagnac + I + T
P_corr = 20,149,603.512 + -132,815.465 - -5.640 + 9.999 + 0.000
P_corr = 20,282,408.839 m
residual = 20,282,408.839 - 20,282,408.337

Modeled/corrected

20,282,408.839 m

Measured raw

20,149,603.512 m

Pre-fit residual

0.502 m

Interpretation: Residual is expected before estimation. The clock term is the largest contributor here, atmospheric terms are meter-level, and the corrected value feeds the linearized model.

Core corrections

Transmit-time estimation and corrected range

From raw observations to a corrected geometric range. Each correction is computed numerically so you can see how meters are gained or lost before solving for position.

Transmit time

ttx≈trx−P/ct_{tx} \approx t_{rx}-P/c
t_rx378,000.000 s
P/c0.067211843 s
t_tx377,999.932345133 s
Interpretation: The first transmit-time estimate moves the signal from reception back to satellite emission.

Satellite clock correction

δts=af0+af1(ttx−Toc)+af2(ttx−Toc)2\delta t_s=a_{f0}+a_{f1}(t_{tx}-T_{oc})+a_{f2}(t_{tx}-T_{oc})^2
delta t_s4.430247e-4 s
c delta t_s132,815.465 m
range sign used-132,815.465 m
Interpretation: The exported correction is applied with the observation-model sign convention.

Broadcast ephemeris snippet

RINEX navigation record

G09 2024 11 14 10 00 00  4.430247051307D-04
     toe = 378000.000 s       sqrtA/e/i0 exported in NAV
     t_tx = 377999.932345 s   delta_t_s = 4.430247051307D-04 s
     Sagnac = -5.639826 m     I_L1 = 9.999084 m
Earth rotation Sagnac correction geometry

Earth rotation / Sagnac

Δρrot=ωEc(xsyr−ysxr)\Delta\rho_{rot}=\frac{\omega_E}{c}(x_s y_r-y_s x_r)
x_s, y_s8,585,785.9, 17,860,591.7
x_r, y_r3,246,374.6, 4,052,665.4
Delta rho rot-5.640 m

Ionosphere correction

Klobuchar model correction

modelKlobuchar
elevation74.43 deg
I9.999 m
Interpretation: The L1 code ionosphere correction is meter-level for this satellite.

Troposphere correction

Troposphere model correction

modelExported correction
elevation74.43 deg
T0.000 m
Interpretation: This Qelora export reports zero troposphere for the selected processing configuration.

Numerical checkpoint

Corrected geometric range

ρcorr=P+cδts−Δρrot−I−T\rho_{corr}=P+c\delta t_s-\Delta\rho_{rot}-I-T
rho_corr = 20,149,603.512 + -132,815.465 - (-5.640) + 9.999 + 0.000
rho_corr = 20,282,408.839 m

Summary of corrections

Linearized measurement model

Turn corrected ranges into a solvable system

δρ=Aδx+b+ϵ\delta\rho=A\delta x+b+\epsilon

State correction vector

δx=[δxδyδzδtr]T\delta x=\begin{bmatrix}\delta x&\delta y&\delta z&\delta t_r\end{bmatrix}^{T}

Design matrix

As=[−ex−ey−ez1]A_s=\begin{bmatrix}-e_x&-e_y&-e_z&1\end{bmatrix}

Bias/correction vector

b=[b1b2⋯bm]Tb=\begin{bmatrix}b_1&b_2&\cdots&b_m\end{bmatrix}^{T}

delta rho: observation residuals. b: modeled corrections. epsilon: measurement noise.

Geometry insight

Design matrix row and geometric residual (G09)

Receiver to satellite line-of-sight geometry

Design matrix row

e[-0.263254, -0.680782, -0.683544]
A_s[-0.263254, -0.680782, -0.683544, 1.000000]
A_clock1.000000

Geometric range

Receiver approximate[3,246,374.6, 4,052,665.4, 3,693,176.1]
Satellite ECEF[8,585,785.9, 17,860,591.7, 17,557,100.6]
rho_geom = ||s-r||20,282,408.337 m
Corrected observation20,282,408.839 m
Residual0.502 m

Geometric residual

0.5021 m

Build all satellite rows

Every satellite becomes one row of H

SVElevPseudorange correctedGeometric rangeResidualDesign matrix row A_s
G0974.4320,282,408.83920,282,408.3370.5021[-0.2633, -0.6808, -0.6835, 1.0000]
G0765.6021,128,343.94421,128,344.609-0.6655[-0.6290, -0.2606, -0.7324, 1.0000]
G0447.1021,469,331.00521,469,331.378-0.3731[0.1118, -0.9352, -0.3359, 1.0000]
G3036.9022,434,075.29222,434,075.1980.0939[-0.9372, 0.1088, -0.3315, 1.0000]
G1619.6023,425,789.90123,425,789.970-0.0691[0.5779, -0.2982, -0.7597, 1.0000]

Pre-fit residuals

m

5

mean

-0.102 m

RMS

0.412 m

Weighted least squares

Solve the linearized system

Δθ=(HTWH)−1HTWΔρ\Delta\theta=(H^TWH)^{-1}H^TW\Delta\rho
wi=(sin⁡(Ei)/σ0)2,σ0=0.3 mw_i=(\sin(E_i)/\sigma_0)^2,\quad \sigma_0=0.3\ m
Hdesign matrix
Wweight matrix
Delta rhopre-fit residual vector
Delta thetastate update

Solution (this iteration)

Iteration 1

Delta x

-17.9810 m

Delta y

17.9596 m

Delta z

5.5403 m

Delta cdt

-0.1023 m

Position update norm

26.0108 m

Satellites used

5

Interpretation: All exported measurements are used. Model: pseudorange SPP. Weights: elevation-based in the data schema.

Solve and iterate

Repeat until the update collapses

1

Form linearized model

Build H and Delta rho

2

Weighted least squares

Compute Delta theta

3

Update state

xk+1=xk+Δxx_{k+1}=x_k+\Delta x
4

Check convergence

Position & clock norms

5

Repeat until solved

Thresholds met

Final iteration

2

Position update

0.00009981 m

RMS post-fit

0.398870 m

Converged

Yes

Quality checks

Residual diagnostics, DOP, and matrix structure

Iteration history

itupdateRMSclock
126.0107650.398873-0.102312
20.0001000.398870-0.102321

Convergence chart

Converged at iteration 2.

Post-fit residuals

All residuals are within +/-1 m in this exported solution.

DOP summary

GDOP

3.69

PDOP

3.15

HDOP

2.12

VDOP

2.34

TDOP

1.93

Design matrix heatmap

A_xA_yA_zA_clockG09-0.26-0.68-0.681.00G07-0.63-0.26-0.731.00G040.11-0.94-0.341.00G30-0.940.11-0.331.00G160.58-0.30-0.761.00

Prefit to postfit check

Worked example summary

Final SPP solution

ECEF position

X

3,246,356.6055 m

Y

4,052,683.3535 m

Z

3,693,181.6012 m

Geodetic

Latitude

35.60388432 deg

Longitude

51.30391819 deg

Height

1,092.738 m

Metrics

RMS post-fit

0.398870 m

PDOP

3.153

Used satellites

5

Model

code

SPP solution complete

You built this

SPP solution complete

Receiver position from pseudorange observations

You successfully computed a GNSS receiver position from pseudorange observations using Single Point Positioning. Next up: examine DOP and error sources, then build a carrier-phase relative-positioning model with double differences.

Receiver final position illustration

Continue your learning

Next numerical threads

DOP and Geometry

Understand GDOP, PDOP and satellite geometry

Open

Error Sources in GNSS

Ionosphere, troposphere, multipath and more

Open

LAMBDA Ambiguity Resolution

Search integer ambiguities using covariance and decorrelation

Open

Double Difference Positioning

Build a correlated carrier-phase model from real RINEX data

Open

Real Data Walkthroughs

RINEX examples and case studies

Open

Questions this example answers

Why are at least four satellites normally needed for SPP?

A three-dimensional code solution has four basic unknowns: X, Y, Z and receiver clock bias. Four independent measurements are a minimum, not a guarantee of a well-conditioned solution. Additional satellites provide redundancy and support residual checks.

Why calculate satellite position at transmission time?

The received signal left the moving satellite earlier. Use an estimated travel time and clock information to evaluate the satellite state at transmission, then handle Earth rotation consistently when comparing it with reception-frame receiver coordinates.

Does a small least-squares residual prove an accurate position?

No. A residual measures fit to the chosen observation model. Unmodelled biases, weak geometry and stochastic assumptions can leave a small residual without proving agreement with independently known coordinates.

Primary sources and further reading

The specifications and research below support the model and terminology. Numerical exports and teaching assumptions are documented in the lesson itself.

Browse the scientific reference library →

Continue with the numbers

GNSS terms and units · Report a correction