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

The three exits of a gated geodetic iteration The iteration starts from an initial estimate and computes a correction. Three exits follow. If the correction is below the convergence threshold the result is returned. If the iteration count reaches its cap the routine raises rather than returning the last estimate, because an unconverged coordinate that looks converged is the failure mode this gate exists to prevent. If the correction grows instead of shrinking, the geometry is near-antipodal and the routine says so explicitly. Initial estimate Compute correction one iteration |correction| < eps ? iteration == max ? Return result converged Raise unconverged or near-antipodal no yes yes no

Figure — a gated iteration has three exits, and only one of them returns a coordinate.

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.

Size of the correction at each iteration Bar chart of the remaining correction in metres at each of five iterations: 0.24, 0.00031, 5.2e-10, and then at the floor of double precision. The correction is roughly squared each pass, which is what quadratic convergence looks like, so a cap of five iterations is generous for well-conditioned geometry and a cap that is actually reached signals a geometry problem, not a tuning problem. 1e-12 1e-11 1e-10 1e-09 1e-08 1e-07 1e-06 1e-05 0.0001 0.001 0.01 0.1 1 m 2.4e-01 iter 1 3.1e-04 iter 2 5.2e-10 iter 3 1.0e-12 iter 4 1.0e-12 iter 5

Figure — quadratic convergence: the correction squares each pass, so five iterations is a generous cap.

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.

Frequently Asked Questions

Why does the iteration fail near antipodal points?

Because the geometry becomes genuinely ambiguous there. For two nearly antipodal points on an ellipsoid, many geodesics of almost equal length connect them, and the iterative formulation oscillates between them instead of converging. This is a property of the problem, not a bug in the implementation, and it is why the iteration must be capped and must raise rather than return its last estimate. Modern series solutions handle the case directly; the classical iteration does not.

What convergence threshold should I use?

Tie it to the ground precision you actually need rather than to a round number. A threshold of 1e-12 radians in latitude corresponds to roughly six micrometres on the ground, comfortably below any survey tolerance and still reached in three or four iterations for well-conditioned geometry. Pushing to the last representable bit costs iterations, buys nothing measurable and can prevent convergence altogether when the correction bottoms out in floating-point noise.

Should I return the number of iterations to the caller?

Yes. It costs nothing and it is the cheapest diagnostic you will ever have: a routine that normally converges in three iterations and suddenly takes nine is telling you something about the data it was given, and that signal is invisible if only the coordinate is returned. Log it with the result in the audit record.