CoordTransform Cadastral & High-Precision Python

Python Script for Geodetic Inverse Problem Solving: Deterministic Vincenty with Antipodal Fallback

Solving the geodetic inverse problem — the ellipsoidal distance, forward azimuth, and reverse azimuth between two known coordinates — must return a finite, auditable result within ±0.001 m of the reference solution even at antipodal configurations where the classical iteration refuses to converge. This operation feeds directly into error distribution modeling in Python, the parent topic under Algorithmic Math & Geodetic Workflows, because every distance and bearing it produces becomes an observation whose residual is later propagated through a covariance matrix.

Deterministic Vincenty inverse-solve decision logic A flowchart of the geodetic inverse solver: coincident points short-circuit to a zero result; otherwise the Vincenty iteration refines lambda on the auxiliary sphere and tests the change against a 1e-12 radian tolerance. Convergence routes to the ellipsoidal delta-sigma series; reaching the 1000-iteration cap routes to a spherical great-circle fallback that always returns a finite distance with a flag set. Geodetic inverse solve P1(φ₁,λ₁) → P2(φ₂,λ₂) Coincident? Yes return d = 0, az = 0 No Reduced lat U₁, U₂; λ = L Vincenty step: recompute σ, update λ on auxiliary sphere |Δλ| < 1e-12 rad? No Yes: converged iter ≥ 1000? No — re-iterate Yes Great-circle fallback → finite dist, flag = True Ellipsoidal Δσ series → distance, fwd / rev azimuth

Why the Iteration Must Be Gated

The Vincenty inverse formula evaluates the angular separation σ\sigma between two reduced-latitude points and refines the longitude difference λ\lambda on the auxiliary sphere until it stops changing. Each pass computes

sin⁡σ=(cos⁡U2sin⁡λ)2+(cos⁡U1sin⁡U2−sin⁡U1cos⁡U2cos⁡λ)2 \sin\sigma = \sqrt{(\cos U_2 \sin\lambda)^2 + (\cos U_1 \sin U_2 - \sin U_1 \cos U_2 \cos\lambda)^2}

and updates λ\lambda from the previous iterate. The loop is declared converged when ∣λi−λi−1∣<10−12|\lambda_i - \lambda_{i-1}| < 10^{-12} radians, a threshold tight enough to hold sub-millimetre distance stability on the WGS84 ellipsoid. For nearly antipodal points the update oscillates instead of contracting, so an unbounded loop will spin until it emits NaN or hangs the pipeline. A hard iteration cap converts that pathology into a deterministic branch: when the cap is reached, the solver routes to a spherical great-circle approximation that always returns a finite value, flagging the result so downstream stages can widen the uncertainty buffer. This is the same numerical-stability discipline that governs least-squares adjustment for control networks, where a single non-finite observation invalidates the entire normal-equation system.

A second hazard is IEEE-754 micro-drift: cos⁡2α=1−sin⁡2α\cos^2\alpha = 1 - \sin^2\alpha can dip fractionally below zero from rounding, and the subsequent math.sqrt or division then raises a domain error. Clamping cos⁡2α\cos^2\alpha to a non-negative floor before it propagates keeps the post-convergence ellipsoidal terms well-defined.

Complete Runnable Implementation

The solver below is self-contained, type-hinted, and frozen-dataclass output for safe reuse inside automated pipelines. It hard-caps iteration, clamps the eccentricity term, and falls back to a great-circle solution when the ellipsoidal series fails to contract.

from __future__ import annotations
import math
from dataclasses import dataclass


@dataclass(frozen=True)
class GeodeticInverseResult:
    distance_m: float
    forward_azimuth_deg: float
    reverse_azimuth_deg: float
    convergence_iterations: int
    fallback_triggered: bool


class GeodeticInverseSolver:
    """Deterministic Vincenty inverse solver with antipodal fallback routing.

    Implements the NGS Vincenty inverse formulation on the WGS84 ellipsoid
    with an explicit iteration cap and a survey-grade convergence tolerance.
    """

    # WGS84 ellipsoid parameters (EPSG:7030)
    A: float = 6378137.0                # semi-major axis (m)
    F: float = 1.0 / 298.257223563      # flattening
    B: float = A * (1.0 - F)            # semi-minor axis (m)

    MAX_ITERATIONS: int = 1000
    CONVERGENCE_TOLERANCE: float = 1e-12  # radians on lambda

    @classmethod
    def solve(cls, lat1_deg: float, lon1_deg: float,
              lat2_deg: float, lon2_deg: float) -> GeodeticInverseResult:
        # Short-circuit coincident coordinates before the iteration runs.
        if (math.isclose(lat1_deg, lat2_deg, abs_tol=1e-12)
                and math.isclose(lon1_deg, lon2_deg, abs_tol=1e-12)):
            return GeodeticInverseResult(0.0, 0.0, 0.0, 0, False)

        lat1, lon1, lat2, lon2 = map(math.radians,
                                     (lat1_deg, lon1_deg, lat2_deg, lon2_deg))
        U1 = math.atan((1.0 - cls.F) * math.tan(lat1))  # reduced latitude
        U2 = math.atan((1.0 - cls.F) * math.tan(lat2))
        L = lon2 - lon1
        lambda_val = L

        sinU1, cosU1 = math.sin(U1), math.cos(U1)
        sinU2, cosU2 = math.sin(U2), math.cos(U2)

        fallback_triggered = False
        iterations = 0
        sin_sigma = cos_sigma = sigma = cos2_sigma_m = cos_sq_alpha = 0.0
        sin_lambda = cos_lambda = 0.0

        for _ in range(cls.MAX_ITERATIONS):
            iterations += 1
            sin_lambda = math.sin(lambda_val)
            cos_lambda = math.cos(lambda_val)

            sin_sigma = math.sqrt(
                (cosU2 * sin_lambda) ** 2
                + (cosU1 * sinU2 - sinU1 * cosU2 * cos_lambda) ** 2
            )
            if sin_sigma == 0.0:  # coincident on the auxiliary sphere
                return GeodeticInverseResult(0.0, 0.0, 0.0, iterations, False)

            cos_sigma = sinU1 * sinU2 + cosU1 * cosU2 * cos_lambda
            sigma = math.atan2(sin_sigma, cos_sigma)
            sin_alpha = (cosU1 * cosU2 * sin_lambda) / sin_sigma
            cos_sq_alpha = max(1.0 - sin_alpha ** 2, 0.0)  # clamp IEEE-754 drift

            cos2_sigma_m = (cos_sigma - (2.0 * sinU1 * sinU2) / cos_sq_alpha
                            if cos_sq_alpha != 0.0 else 0.0)
            C = (cls.F / 16.0) * cos_sq_alpha * (4.0 + cls.F * (4.0 - 3.0 * cos_sq_alpha))

            lambda_prev = lambda_val
            lambda_val = L + (1.0 - C) * cls.F * sin_alpha * (
                sigma + C * sin_sigma * (
                    cos2_sigma_m + C * cos_sigma * (-1.0 + 2.0 * cos2_sigma_m ** 2)
                )
            )
            if abs(lambda_val - lambda_prev) < cls.CONVERGENCE_TOLERANCE:
                break
        else:
            fallback_triggered = True

        if fallback_triggered:
            # Spherical great-circle fallback: always finite near antipodal points.
            cos_sig = math.sin(lat1) * math.sin(lat2) + \
                math.cos(lat1) * math.cos(lat2) * math.cos(L)
            sigma = math.acos(max(-1.0, min(1.0, cos_sig)))
            distance_m = cls.A * sigma
            fwd = math.atan2(math.sin(L) * math.cos(lat2),
                             math.cos(lat1) * math.sin(lat2)
                             - math.sin(lat1) * math.cos(lat2) * math.cos(L))
            rev = math.atan2(math.sin(-L) * math.cos(lat1),
                             math.cos(lat2) * math.sin(lat1)
                             - math.sin(lat2) * math.cos(lat1) * math.cos(-L))
            return GeodeticInverseResult(
                distance_m=distance_m,
                forward_azimuth_deg=math.degrees(fwd) % 360.0,
                reverse_azimuth_deg=(math.degrees(rev) + 180.0) % 360.0,
                convergence_iterations=iterations,
                fallback_triggered=True,
            )

        # Post-convergence ellipsoidal distance via the Vincenty series.
        u_sq = cos_sq_alpha * (cls.A ** 2 - cls.B ** 2) / cls.B ** 2
        A_c = 1.0 + (u_sq / 16384.0) * (4096.0 + u_sq * (-768.0 + u_sq * (320.0 - 175.0 * u_sq)))
        B_c = (u_sq / 1024.0) * (256.0 + u_sq * (-128.0 + u_sq * (74.0 - 47.0 * u_sq)))
        delta_sigma = B_c * sin_sigma * (
            cos2_sigma_m + (B_c / 4.0) * (
                cos_sigma * (-1.0 + 2.0 * cos2_sigma_m ** 2)
                - (B_c / 6.0) * cos2_sigma_m * (-3.0 + 4.0 * sin_sigma ** 2)
                * (-3.0 + 4.0 * cos2_sigma_m ** 2)
            )
        )

        distance_m = cls.B * A_c * (sigma - delta_sigma)
        fwd = math.atan2(cosU2 * sin_lambda, cosU1 * sinU2 - sinU1 * cosU2 * cos_lambda)
        rev = math.atan2(cosU1 * sin_lambda, -sinU1 * cosU2 + cosU1 * sinU2 * cos_lambda)

        return GeodeticInverseResult(
            distance_m=distance_m,
            forward_azimuth_deg=math.degrees(fwd) % 360.0,
            reverse_azimuth_deg=(math.degrees(rev) + 180.0) % 360.0,
            convergence_iterations=iterations,
            fallback_triggered=False,
        )

Parameter Reference

Name Type Units Valid range Notes
lat1_deg, lat2_deg float degrees −90 to +90 Geodetic latitudes of the two stations
lon1_deg, lon2_deg float degrees −180 to +180 Geodetic longitudes
CONVERGENCE_TOLERANCE float radians 1e-12 typical Δλ contraction threshold per iteration
MAX_ITERATIONS int count 1000 default Hard cap that triggers the fallback branch
distance_m float metres ≥ 0 Ellipsoidal (or great-circle, on fallback) distance
forward_azimuth_deg float degrees 0 ≤ x < 360 Bearing at point 1 toward point 2
reverse_azimuth_deg float degrees 0 ≤ x < 360 Bearing at point 2 toward point 1
fallback_triggered bool — — True when the great-circle path was used

Worked Example

Solving the inverse between two U.S. National Spatial Reference System control points — Flinders Peak and Buninyong in the classic Vincenty test pair — returns the published ellipsoidal distance to sub-millimetre agreement.

# Flinders Peak (-37°57'03.72030", 144°25'29.52440") to
# Buninyong   (-37°39'10.15610", 143°55'35.38390") in decimal degrees.
r = GeodeticInverseSolver.solve(-37.95103342, 144.42486789,
                                -37.65282114, 143.92649553)
print(round(r.distance_m, 3), round(r.forward_azimuth_deg, 5), r.fallback_triggered)
# 54972.271 306.86816 False

Validation Check

Gate the result against survey tolerance before it is committed to a control network adjustment or a cadastral database.

EXPECTED_M = 54972.271
assert math.isfinite(r.distance_m), "non-finite distance"
assert 0.0 <= r.forward_azimuth_deg < 360.0, "azimuth out of bounds"
assert abs(r.distance_m - EXPECTED_M) < 0.001, "distance drift exceeds survey tolerance"

Common Mistakes

  • Treating a fallback result as ellipsoidal. When fallback_triggered is True, the distance is a spherical great-circle value that can deviate by tens of metres from the true geodesic. Always inspect the flag and inflate the observation variance — or switch to Karney’s algorithm — before passing the value into a covariance propagation step.
  • Omitting the iteration cap. A while abs(dlambda) > tol loop with no bound will hang or emit NaN on antipodal pairs. The for ... else construction here guarantees a deterministic exit on every input.
  • Forgetting the reduced latitude. Vincenty operates on the reduced (parametric) latitude U=arctan⁡((1−f)tan⁡ϕ)U = \arctan((1-f)\tan\phi), not the geodetic latitude. Passing ϕ\phi directly produces a systematic bias that grows with the flattening term — the same class of error that distorts an affine transformation fit for a local grid when the input frame is mislabelled.