Taking altitude into account when calculating geod

Question

i´m currently dealing with gps data combined with precise altitude measurement. I want to calculate the distance between two consecuting points. There is a lot of information out there about calculating distance between two points using the WGS84 ellipsoid and so on.

however, i did not find any information that takes Altitude changes into account for this distance calculation.

does anyone know about some websites, papers, books etc. that describes such a method? thanks

edit: Sql Server 2008 geographic extensions also neglect altitude information when calculating distance.

Answer 1

For starters, you need a model that tells you how the altitude changes on the line between the two points. Without such a model, you don't have any consistent definition of the distance between two points.

If you had a linear model (traveling 50% of the distance between the points also means you went upwards through 50% of the altitude), then you can probably pretend that the entire thing was a right-triangle; i.e. you act as though the world is flat for purposes of determining how the altitude shift affects the distance. The distance along the ground is the base, the altitude change is the height of the triangle, and the hypotenuse is your estimated true travel distance from point to point.

If you want to refine that further, then you can note that the model above is perfectly good for infinitesimal distances, which means that you can iterate across individual deltas of the distance, calculus-style, each time using the current altitude to compute the ground distance and then using the same trigonometric ratio to compute the altitudinal-change contribution to the distance traveled. I'd probably do this in a for() loop with 10 to 100 pieces of the segment, and possibly by trial and error figure out the number of pieces required to get within epsilon of the true value. It would also be possible to work out the line integral to figure out the actual distance between the two points under this model.

Answer 2

You likely don't care about altitude for large 2D distance separatiions. So if the dist you get is over say 20 (or perhaps 50)km, then who cares about the altitude diff (depends on your needs case). Under say 20km, feed in the simple pythagorean addition to the altitude difference. Feed it in smoothly.

Distance between two geo-points?

Answer 3

I implemented a WGS84 distance function using the average of the start and end altitude as the constant altitude. If you are certain that there will be relatively little altitude variation along your path this works acceptably well (error is relative to the altitude difference of your two LLA points).

Here's my code (C#):

    /// <summary>
    /// Gets the geodesic distance between two pathpoints in the current mode's coordinate system
    /// </summary>
    /// <param name="point1">First point</param>
    /// <param name="point2">Second point</param>
    /// <param name="mode">Coordinate mode that both points are in</param>
    /// <returns>Distance between the two points in the current coordinate mode</returns>
    public static double GetGeodesicDistance(PathPoint point1, PathPoint point2, CoordMode mode) {
        // calculate proper geodesics for LLA paths
        if (mode == CoordMode.LLA) {
            // meeus approximation
            double f = (point1.Y + point2.Y) / 2 * LatLonAltTransformer.DEGTORAD;
            double g = (point1.Y - point2.Y) / 2 * LatLonAltTransformer.DEGTORAD;
            double l = (point1.X - point2.X) / 2 * LatLonAltTransformer.DEGTORAD;

            double sinG = Math.Sin(g);
            double sinL = Math.Sin(l);
            double sinF = Math.Sin(f);

            double s, c, w, r, d, h1, h2;
            // not perfect but use the average altitude
            double a = (LatLonAltTransformer.A + point1.Z + LatLonAltTransformer.A + point2.Z) / 2.0;

            sinG *= sinG;
            sinL *= sinL;
            sinF *= sinF;

            s = sinG * (1 - sinL) + (1 - sinF) * sinL;
            c = (1 - sinG) * (1 - sinL) + sinF * sinL;

            w = Math.Atan(Math.Sqrt(s / c));
            r = Math.Sqrt(s * c) / w;
            d = 2 * w * a;
            h1 = (3 * r - 1) / 2 / c;
            h2 = (3 * r + 1) / 2 / s;

            return d * (1 + (1 / LatLonAltTransformer.RF) * (h1 * sinF * (1 - sinG) - h2 * (1 - sinF) * sinG));
        }

        PathPoint diff = new PathPoint(point2.X - point1.X, point2.Y - point1.Y, point2.Z - point1.Z, 0);
        return Math.Sqrt(diff.X * diff.X + diff.Y * diff.Y + diff.Z * diff.Z);
    }

In practice we've found that the altitude difference rarely makes a large difference, our paths are typically 1-2km long with altitude varying on the order of 100m and we see about ~5m change on average versus using the WGS84 ellipsoid unmodified.

Edit:

To add to this, if you do expect large altitude changes, you can convert your WGS84 coordinates to ECEF (earth centered earth fixed) and evaluate straight-line paths as shown at the bottom of my function. Converting a point to ECEF is simple to do:

    /// <summary>
    /// Converts a point in the format (Lon, Lat, Alt) to ECEF
    /// </summary>
    /// <param name="point">Point as (Lon, Lat, Alt)</param>
    /// <returns>Point in ECEF</returns>
    public static PathPoint WGS84ToECEF(PathPoint point) {
        PathPoint outPoint = new PathPoint(0);

        double lat = point.Y * DEGTORAD;
        double lon = point.X * DEGTORAD;
        double e2 = 1.0 / RF * (2.0 - 1.0 / RF);
        double sinLat = Math.Sin(lat), cosLat = Math.Cos(lat);

        double chi = A / Math.Sqrt(1 - e2 * sinLat * sinLat);
        outPoint.X = (chi + point.Z) * cosLat * Math.Cos(lon);
        outPoint.Y = (chi + point.Z) * cosLat * Math.Sin(lon);
        outPoint.Z = (chi * (1 - e2) + point.Z) * sinLat;

        return outPoint;
    }

Edit 2:

I was asked about some of the other variables in my code:

// RF is the eccentricity of the WGS84 ellipsoid
public const double RF = 298.257223563;

// A is the radius of the earth in meters
public const double A = 6378137.0;

LatLonAltTransformer is a class I used to convert from LatLonAlt coordinates to ECEF coordinates and defines the constants above.

Answer 4

In order to do this the first issue you have to address is how to define change in altitude. The normal equations work because they are on a two dimensional surface, however adding the third dimension means that the simple definition of shortest distance is no longer applicable, for example now the thrid dimension is 'in play' your shortest distance could cut through the original ellipsoid. It's a bit quick and dirty, but your best solution might be to assume that the rate of change of alltitude is constant along the original 2D path on the ellipsoid. You can then calculate the 2D distance as a length, work out the rate of change of altitude and then simply use Pythagoras to calculate the increase in length with one side of the triangle being the 2D distance, and the altitude being the second length.

Answer 5

I would suggest that over any distance where using the WGS84 would give you significantly better accuracy that the difference in altitude won't matter. And over any distance where the difference in altitudes matters you should probably just use straight line approximation.