Geodesics on the ellipsoid, and why they are hard

A great circle is the intersection of a sphere with a plane through its centre, which makes it simple to generate, simple to measure and simple to reason about.

On an ellipsoid none of that survives. The shortest path is not a plane curve, its length is an elliptic integral, and it does not return to its starting point after one circuit. It is the same question — what is the shortest route — and it becomes a genuinely harder object.

What survives #

One thing, and it is what makes the problem tractable at all.

Clairaut’s relation holds on any surface of revolution, so it holds here — but in terms of the parametric latitude $\beta$, where $\tan\beta = (1-f)\tan\varphi$:

$$\cos\beta\,\sin\alpha = \text{constant}$$

That is the only conserved quantity, and every practical method for ellipsoidal geodesics is built on it. Vincenty’s begins by converting both endpoints to parametric latitude precisely so that the relation takes its spherical form, then applies a series correction in the eccentricity.

So the auxiliary latitudes are load-bearing: the parametric one is the change of variable that makes an ellipsoid’s geodesics look spherical, and choosing a different member of the family does not make the calculation slightly worse — it makes it answer a different question.

What does not survive #

Three properties that make a great circle easy, all gone.

It is not a plane curve. A geodesic on an ellipsoid does not lie in any plane. It winds, and the plane through its endpoints and the centre is not it.

It does not close. Follow a great circle round a sphere and it returns to where it started. Follow an ellipsoidal geodesic all the way round and it comes back displaced in longitude — the amount depends on the flattening and on the inclination, and the path fills a band rather than tracing a curve.

Its length is an elliptic integral. There is no elementary formula. Every method is either a series in the flattening or an iteration, and the choice between them is about which failure mode is acceptable.

Two methods, and why both are needed #

Vincenty’s inverse method, from 1975, is the classical answer: an iteration on one variable that converges in a handful of steps for ordinary journeys, with a series correction for the distance.

The problem with an iteration is that its accuracy is not visible. It converges to something, and if the series coefficients were wrong it would converge to something wrong with equal confidence. Comparing it against a longer version of the same series would restate rather than test.

So the site computes the geodesic a second way, sharing nothing with the first: it integrates the geodetic reckoning equations forward from the start point with a fourth-order Runge–Kutta step,

$$\frac{\mathrm{d}\varphi}{\mathrm{d}s} = \frac{\cos\alpha}{M(\varphi)}, \qquad \frac{\mathrm{d}\lambda}{\mathrm{d}s} = \frac{\sin\alpha}{N(\varphi)\cos\varphi}, \qquad \frac{\mathrm{d}\alpha}{\mathrm{d}s} = \frac{\sin\alpha\tan\varphi}{N(\varphi)}$$

using the two radii of curvature directly. Vincenty gives a distance and a departure azimuth; integrating from the same point, on that azimuth, for that distance, must arrive at the destination.

Across five journeys of six to eleven thousand kilometres, the worst discrepancy is 29 micrometres.

That is the check worth having. Two independent routes — an iteration on a transformed variable, and a differential equation stepped forward — landing within a fraction of a millimetre of each other over ten thousand kilometres.

How much the flattening is worth #

The numbers the essays need, and they are more interesting than a single correction factor.

Four routes longer, one shorter, and the exception is the informative one.

Cape Town to London runs almost due north. An oblate Earth is flatter pole to pole than a sphere of the mean radius, so a meridional route is shorter than the spherical estimate. Every other route here has a substantial east–west component and runs where the Earth is fatter than the mean sphere, so it is longer.

There is therefore no correction factor. A rule saying “add a quarter of a per cent to spherical distances” is wrong in sign for a third of the journeys on this site, and the sign depends on the bearing rather than on anything about the endpoints.

Why the site uses a sphere anyway #

Stating the justification as a measurement rather than as a convention, since the essays make the assumption everywhere.

The corrections are between 0.20% and 0.36%. The site asserts both bounds: the difference must be under one per cent, or the spherical arguments would need qualifying rather than a footnote; and it must be over 0.01%, or one of the two computations is not using the ellipsoid it claims.

Both halves matter. An assertion that only checked the difference was small would pass if the ellipsoidal routine had quietly reverted to a sphere.

And the structure is identical on both bodies. That the shortest route is not the constant-bearing route, that it runs poleward, that the excess grows with latitude, that Clairaut’s product is conserved — all of those are consequences of rotational symmetry and positive curvature, and none depends on the curvature being constant. The numbers shift by a third of a per cent; the argument does not shift at all.

Which is the correct division of labour. Use a sphere for the shape of the answer, and an ellipsoid when a number is going to be compared with a measurement on the ground.

Where the methods break #

Vincenty’s iteration has a known failure and it is worth naming, because it is the reason the modern algorithm exists.

For nearly antipodal points — two places on opposite sides of the Earth — the iteration converges slowly or not at all. The reason is geometric rather than numerical: near the antipode there are many nearly-shortest paths of almost the same length, the problem is close to degenerate, and an iteration looking for one solution has nothing to hold onto.

On a sphere the antipodal case is degenerate too — every great circle through a point passes through its antipode, so all paths are equally short — but the degeneracy is exact and obvious. On an ellipsoid it is approximate and messy: there is a unique shortest path, and it is surrounded by a family of others within metres of its length.

Karney’s 2013 algorithm solves this by reformulating the problem so that it converges everywhere, including antipodal cases, and to nanometre accuracy. It is what current libraries use, and the site reports whether Vincenty converged rather than assuming it did.

The direct problem, and why it is easier #

Two questions are usually distinguished, and the asymmetry between them says something about the shape of the difficulty.

The direct problem: given a start point, a bearing and a distance, where does the path end? The inverse problem: given two points, what is the distance and the bearing between them?

On a sphere both are a few lines of trigonometry and neither is harder. On an ellipsoid the direct problem is genuinely easier: it is an initial-value problem, the equations can be stepped forward, and the answer is unique and well-conditioned everywhere. The site’s independent check is exactly a direct-problem solver, and it is a fourth-order Runge–Kutta integration of three coupled equations — twenty lines, no iteration, no special cases.

The inverse problem is a boundary-value problem, and boundary-value problems are where the trouble lives. It has to be solved by iteration, the iteration can fail to converge, and near-antipodal endpoints make the solution nearly non-unique.

That asymmetry is why the site can check Vincenty’s inverse method with an integrator at all. The two problems are inverse to each other, so a direct solver run on the inverse solver’s output is a genuine test, and the easier of the two is the one available to do the testing.

What “shortest” stops meaning #

One property of ellipsoidal geodesics deserves stating because it is genuinely surprising and it has no spherical analogue.

A geodesic is locally shortest: any short piece of it is the shortest route between its own endpoints. It is not always globally shortest. Past a certain distance — the point where nearby geodesics from the same origin begin to cross — the path stops being the shortest route even though it is still a geodesic.

On a sphere that happens at the antipode and only there, which is why the distinction rarely comes up. On an ellipsoid the crossing locus is not a point but a curve, and its position depends on the departure bearing.

So “the geodesic between two points” is slightly the wrong phrase. There are infinitely many geodesics between any two points, one of them is the shortest, and identifying it is part of what the inverse problem does. Vincenty’s method returns the shortest one when it converges, and near the antipode the several candidates it is choosing between are within metres of each other — which is the same fact as the convergence failure, described geometrically instead of numerically.

What a geodesic is for #

The practical uses are narrower than the mathematics suggests, and worth listing because they explain why the accuracy is pursued so far.

Boundaries. A great many international and administrative boundaries are defined as geodesics between coordinates — maritime boundaries especially, where the equidistance line between two coasts is constructed geodesically. A boundary is a legal object and a metre matters.

Aviation and shipping. Route planning uses geodesics as the baseline, modified for weather and airspace. Here a kilometre matters and a metre does not.

Geodetic networks. The connection between two survey stations is a geodesic, and the whole business of adjusting a network is reconciling measured lines with computed ones. This is where nanometre accuracy is not absurd: the computation must be far more accurate than the measurement, so that the residuals are the measurement’s error rather than the model’s.

That last case is why Karney’s accuracy is worth having. It is not that anyone can measure a nanometre; it is that a model whose error is comparable to the observations’ contaminates the adjustment, and pushing the model’s error to nothing removes it from the analysis entirely.

What the map has to do with it #

The connection back to projections, since this is a site about them.

A geodesic on the ellipsoid, projected, is a curve on the map — and whether it is straight is a property of the map. No projection renders ellipsoidal geodesics straight in general: the gnomonic projection straightens spherical great circles exactly, and on an ellipsoid it straightens them only to the extent that the flattening is small.

That is a small crack in one of the site’s cleaner results. The gnomonic projection’s uniqueness — every great circle straight, measured at 3×10⁻¹⁶ — is a theorem about surfaces of constant curvature, and an ellipsoid is not one. The projection remains an excellent approximation for route planning and it is no longer exact.

Which is the same pattern as everywhere in this field: a result that is exact on a sphere becomes approximate on an ellipsoid, with a known error, and the honest thing is to say which body the claim is about.

What a sphere of the right radius cannot fix #

An obvious repair suggests itself and does not work, which is worth showing because it is the first thing anybody tries.

If the corrections are between −0.36% and +0.33%, why not choose a sphere radius that removes them? Pick the radius that makes the mean error zero, and quote spherical distances with a corrected constant.

It fails because the corrections have both signs and their sign depends on the bearing. A radius large enough to fix Cape Town to London makes every east–west route worse; a radius fitted to the east–west routes makes the meridional ones worse. There is no single radius that works, because the discrepancy is not a scale error — it is a shape error, and a sphere has no parameter that can represent a shape.

The best a sphere can do is reduce the worst case, which is the same minimax trade as a grid’s scale factor and buys about a factor of two. The mean radius of 6,371.0088 km used throughout this site is the authalic radius — the radius of the sphere with the same surface area as WGS84 — which is the right choice for area comparisons and is not optimised for distances at all.

Stating which mean was taken matters. There are at least four defensible mean radii for the Earth, they differ by a few kilometres, and a distance quoted without saying which is uncertain at the level the whole discussion is about.

Vincenty’s formulae agree with the integrated reference to micrometres over ordinary journeys, and there is one neighbourhood of every point where they agree with nothing at all.

What was computed here #

Distances on the sphere use the haversine formula. Distances on the ellipsoid use Vincenty’s inverse method on WGS84, with the iteration count and convergence reported rather than assumed.

The check is the independent one. Integrating the geodetic reckoning equations forward from each start point, on Vincenty’s departure azimuth, for Vincenty’s distance, must arrive at the destination — and it does, to 29 micrometres over journeys of up to eleven thousand kilometres. The two computations share no algebra: one is a fixed-point iteration on an auxiliary longitude with a series correction, the other is a fourth-order Runge–Kutta step through three coupled differential equations.

Two bounds are asserted on the flattening’s effect. It must be under one per cent, which licenses the spherical arguments elsewhere on this site. And it must be over 0.01%, which is the check that the ellipsoidal computation is genuinely using an ellipsoid — a routine that had silently reverted to a sphere would pass the first test and fail this one.

What the pictures cannot show #

The geodesic’s non-planarity, which is the thing that makes it hard. A curve on a sphere and a curve on an ellipsoid look identical at any drawing scale, and the difference between them — that one lies in a plane and the other does not — is a fact about the third dimension that a projected figure has already discarded.

The figures also cannot show the antipodal failure. It happens for a narrow band of endpoint pairs, none of the five journeys drawn here is in it, and a picture of an iteration not converging is a picture of nothing.

What the iteration costs a query #

The essay establishes that there is no closed form and that Vincenty’s iteration is the practical answer. The applied field meets the consequence from the other side: a nearest-neighbour query over a large candidate set cannot afford one iteration per candidate, so it is answered in the plane instead — and sometimes answered wrongly.

Over Europe, 5 of 400 queries return a different site when distances are compared in the stored plane rather than on the ellipsoid, with the wrong answer up to 24 kilometres further away. Over the whole world the rate rises to 26 of 400 and the penalty to 1,355 kilometres. The repair is the structure this kind of iteration usually forces: filter with the cheap metric to a width the projection’s scale spread bounds, then decide with the exact one.

The same distinction reaches a line rather than a ranking: a stored straight segment from London to Tokyo and the geodesic the iteration would compute are 2,961 kilometres apart at the worst point, which is a disagreement no tolerance on a query would ever have caught.

Who found it, and when #

Clairaut gave the relation in 1743, working on the figure of the Earth for the same reason the Lapland and Peru expeditions were sent.

Bessel worked out the practical solution in 1825, reducing the ellipsoidal problem to a spherical one by way of the parametric latitude — the substitution every method has used since. Legendre and Gauss both contributed to the theory in the same period, and Gauss’s work on it fed directly into the surface theory that produced the Theorema Egregium.

Thaddeus Vincenty published his iterative formulation in 1975, at the US National Geodetic Survey, explicitly designed for the small computers of the period — the paper’s title mentions “nested equations” because the point was to make the calculation fit. Charles Karney’s reformulation of 2013 removed the antipodal failure and reached nanometre accuracy, and it is what the current libraries implement.

Two hundred and seventy years from Clairaut’s relation to an algorithm that works everywhere, for a question with a one-sentence statement.

Where this goes next #

The invariant that makes the problem solvable is the great-circle vertex. The body it is solved on is the Earth is a sphere, and when it is not. And the spherical version of the whole question is the shortest route is not straight.