What method is used by the mapping toolbox to determine geodesic distances between points within various geodetic reference frames?

I have a radar tracking problem wherin I need to calculate the distance between two points.

While pursuing the answer to this problem I have encountered several methods to acomplish this task.

1. The Vincenty formula was a commonly referenced solution to this problem. It has also been made into a few different Matlab libraries. Michael Kleder and themaze are two users that I found having impletmented the same Vincenty method.
2. Charles Karney published on this topic in 2012 claiming to have a method that can take advantage of the computational power of modern computers. Additionally, this new method improves on the Vincenty method in terms of accuracy and ability to find a solution for points that are nearly antipodal. Karney has also posted a library of functions to accomplish this in matlab.
3. There are several built-in functions in the Mapping Toolbox, such as geodetic2aer or distance, which return the measure of distance between the two inputted points.

There appear to be no references for the mapping toolbox functions which detail what methods are used to calculate these distances. Does anybody know what method the mapping toolbox uses to calculate these things? How accurate is the method matlab uses? Should I opt for the user uploaded methods instead of built-in functions?

The documentation for the distance() function in matlab describes a decrease in accuracy as the distance between points increases and the documentation warns of calculation breakdown for antipodal points. This suggests that the Vincenty formulas was used but it isn't directly said. Any insight is helpful here.

Answer 1

[Jack]

Under the hood, the Mapping Toolbox’s distance (and the related geodetic functions) solve the “inverse geodesic” problem on the reference ellipsoid using the classic Vincenty algorithm. When you supply a sphere, it simply falls back to the great‐circle formula.

Key points:

• Ellipsoid mode: Vincenty’s inverse formula (1975)—accurate to a few tenths of a millimeter over most distances but known to fail or stall for nearly antipodal points.
• Sphere mode: Standard spherical law of cosines or haversine.

Because it really is Vincenty’s routine, you will see the documentation’s warning about degrading accuracy at long distances and breakdown at antipodes. Karney’s 2012 “modern” algorithm (available via his MATLAB file exchange code) fixes those edge‐cases and is numerically more robust, but the built-in Mapping Toolbox hasn’t switched over to it.

If you need guaranteed convergence for any two lat/lon pairs, you should plug in Karney’s routines (or another robust inverse‐geodesic solver) instead of relying on distance(). For most everyday use within a few thousand kilometers, the toolbox’s Vincenty‐based distance() is perfectly adequate.

— Follow me so you can message me anytime with future questions. If this helps, please accept the answer and upvote it as well.

[Steve K]

Thanks for the explanation Jack. It confirms my suspisions at least.

Some follow up questions, that arise from this knowledge:

1. Since the distance() function doesn't take in altitude, does this assume that the inverse geodesic that is solved for lies directly on the ellipsoid surface?
2. Conversely, given that I have points above the surface of the ellipsoid (LLA), am I correct to use the geodetic2aer() function to get a slant range between them, and use the elevation angle to get a geodesic distance along an ellipsoid at height H?

An additional concern, which is probably not appropriate for this forum, but if you have insight I would appreciate it. I have some developers on my team which prefer to use python and the pymap3d package. In that package is a function that is defined very similarly to geodetic2aer(), and it returns the same info. The source code of that package doesn't state which method is used to solve that inverse problem. Do you know what is used or how I can go about finding out? I want to compare apples to apples between code packages to eliminate some disagreements between outputs.

[Jack]

No problem! Here are the answers to your follow up q's:

1. Ellipsoid‐only: The built-in distance always solves on the reference ellipsoid surface—it ignores altitude. In other words, it treats both points as if their height were zero.
2. Slant vs. Ground Range: For true 3D (LLA) separation, use geodetic2aer. It converts each (lat, lon, alt) to ECEF, takes the straight-line (slant) distance, and returns Az/El and slant range. If you need the “geodesic along a parallel at height H,” you can compute

[az,el,slant] = geodetic2aer(lat1,lon1,alt1, lat2,lon2,alt2, wgs84);
groundRange = slant .* cosd(el);

That groundRange is the horizontal distance over an ellipsoid offset by your altitudes.

1. Python’s pymap3d

Their geodetic2aer does essentially the same ECEF → local ENU vector math—it doesn’t run a Vincenty or Karney inverse-geodesic under the hood. You’ll see on github it uses trigonometric ECEF conversions and simple Cartesian differences, not an iterative ellipsoidal solver.

Hope this helps!

[Steve K]

Thanks so much for the clarifications.

I have been using the geodetic2aer() function the way you described. I was doing it based on intuition. So it is nice to have some confirmation. I have been computing the flat range the way you mention as well, I was just unsure if it would be akin to a geodesic or just a straight line that might "cut through" the ellipsoid to a degree.

Concerning pymap3d, are you saying that the methods used by Matlab and pymap3d are the same when it comes to converting to ECEF and/or when transforming from one coordinate system to another? Or, is there some iterative way to properly convert from geodetic to ECEF that pymap3d is not doing?

I bring this up because I get different outputs from both functions for the same inputs. Ex. The slant ranges for the same points are roughly 200m different.

[Jack]

Make sure both tools use the same ellipsoid and the same lat–lon order/units:

• Ellipsoid: MATLAB’s geodetic2ecef defaults to WGS84. In pymap3d you must explicitly request WGS84 (otherwise it may default to a sphere).
• Inputs: Both expect lat, lon in degrees (and altitude in metres). If you swap them or feed radians, you’ll see hundreds of metres of error.

Once you match those, your slant ranges will agree to millimetre-level accuracy.

[Steve K]

We had some labeling problems facepalm

Turns out that a slant range was labeled as ground range in one of our tools. They agree now.

I appreciate the help, at least I've learned a good deal about the back end of both packages now.

Are there any common tripping points I need to consider when converting from LLA to ECEF? I know there are a few methods to do the conversions, checking the wiki at least. I also read a paper published by Olson which suggested that the accuracy difference between methods was on the order of 10^-9m, so I do not think worrying about this is really worth the brain power.

[Jack]

You're right that the different methods generally yield very high accuracy (typically not worth the brain power 🙂‍↕️). However, here are a few common things to be on the lookout for:

• Ellipsoid Parameters: Ensuring both systems are using the exact same ellipsoid definition (semi-major axis, flattening) is crucial.
• Units
• Latitude/Longitude Order
• Sign Conventions
• Numerical Precision: Differences in floating-point precision across different platforms or libraries could theoretically lead to extremely small variations in results, though this is unlikely to be the cause of a 200m difference.

In conclusion, given that you've resolved the labeling issue and the Olson paper suggests negligible accuracy differences between LLA to ECEF conversion methods, it's likely that ensuring consistent ellipsoid parameters and input units/order will be the most important factors for you to manage.