Why is there a difference in SqlGeography STDistance?

Question

I have been tracking down a bug in my code, and I've found that it's because the Microsoft c# SqlGeography 2014 library returns a slightly different result for STDistance than my regular code for calculating the distance between points.

I wrote a small console exe to demonstrate the problem, but I can't figure out why I'm getting such a different result?

static void Main(string[] args) {
        double earthRadius = 6378137; // meters => from both nad83 & wgs84
        var a = new { lat = 43.68151632, lng = -79.61162263 };
        var b = new { lat = 43.67575602, lng = -79.59586143 };

        // sql geography lib
        SqlGeographyBuilder sgb;
        sgb = new SqlGeographyBuilder();
        sgb.SetSrid(4326);
        sgb.BeginGeography(OpenGisGeographyType.Point);
        sgb.BeginFigure(a.lat, a.lng);
        sgb.EndFigure();
        sgb.EndGeography();
        SqlGeography geoA = sgb.ConstructedGeography;

        sgb = new SqlGeographyBuilder();
        sgb.SetSrid(4326);
        sgb.BeginGeography(OpenGisGeographyType.Point);
        sgb.BeginFigure(b.lat, b.lng);
        sgb.EndFigure();
        sgb.EndGeography();
        SqlGeography geoB = sgb.ConstructedGeography;

        // distance cast from SqlDouble
        double geoDistance = (double)geoA.STDistance(geoB);

        // math!
        double d2r = Math.PI / 180; // for converting degrees to radians
        double lat1 = a.lat * d2r,
            lat2 = b.lat * d2r,
            lng1 = a.lng * d2r,
            lng2 = b.lng * d2r,
            dLat = lat2 - lat1,
            dLng = lng2 - lng1,
            sin_dLat_half = Math.Pow(Math.Sin(dLat / 2), 2),
            sin_dLng_half = Math.Pow(Math.Sin(dLng / 2), 2),
            distance = sin_dLat_half + Math.Cos(lat1) * Math.Cos(lat2) * sin_dLng_half;

        // math distance
        double mathDistance = (2 * Math.Atan2(Math.Sqrt(distance), Math.Sqrt(1 - distance))) * earthRadius;

        // haversine
        double sLat1 = Math.Sin(a.lat * d2r),
                sLat2 = Math.Sin(b.lat * d2r),
                cLat1 = Math.Cos(a.lat * d2r),
                cLat2 = Math.Cos(b.lat * d2r),
                cLon = Math.Cos((a.lng * d2r) - (b.lng * d2r)),
                cosD = sLat1 * sLat2 + cLat1 * cLat2 * cLon,
                d = Math.Acos(cosD);

        // math distance
        double methDistance = d * earthRadius;


        // write the outputs
        Console.WriteLine("geo distance:\t" + geoDistance);    // 1422.99560435875
        Console.WriteLine("math distance:\t" + mathDistance);  // 1421.73656776243
        Console.WriteLine("meth distance:\t" + methDistance);  // 1421.73656680185
        Console.WriteLine("geo vs math:\t" + (geoDistance - mathDistance));     // 1.25903659632445
        Console.WriteLine("haversine vs math:\t" + (methDistance - methDistance)); // ~0.00000096058011
    }

Is Microsoft using a different calculation method? Being off by over 1 meter when calculating distances less than 1.5Km is a huge discrepancy.

Answer 1

Ok, so after much digging I found the answer, and Microsoft is more correct.

Specifically, they are using Vincenty's formulae. Accuracy is within 0.5mm (not metre, half a millimetre) instead of 0.3% with Haversine formula.

The reason is that Haversine (used by me, Google, and Bing Maps too apparently) is fast, but relies on a spherical Earth instead of an ellipsoid. Microsoft uses the ellipsoid Earth to calculate distances instead of a sphere providing more accurate results.

I implemented Vincenty's method in c# like this and it's worked so far, but is no where near production ready.

    const double d2r = Math.PI / 180;  // degrees to radians
    const double EARTH_RADIUS = 6378137;  // meters
    const double EARTH_ELLIPSOID = 298.257223563; // wgs84
    const double EARTH_BESSEL = 1 / EARTH_ELLIPSOID;
    const double EARTH_RADIUS_MINOR = EARTH_RADIUS - (EARTH_RADIUS * EARTH_BESSEL); // 6356752.3142 meters => wgs84

    static double vincentyDistance(double lat1, double lng1, double lat2, double lng2) {
        double L = (lng2 - lng1) * d2r,
                U1 = Math.Atan((1 - EARTH_BESSEL) * Math.Tan(lat1 * d2r)),
                U2 = Math.Atan((1 - EARTH_BESSEL) * Math.Tan(lat2 * d2r)),
                sinU1 = Math.Sin(U1),
                cosU1 = Math.Cos(U1),
                sinU2 = Math.Sin(U2),
                cosU2 = Math.Cos(U2),
                lambda = L,
                lambdaP,
                iterLimit = 100,
                sinLambda,
                cosLambda,
                sinSigma,
                cosSigma,
                sigma,
                sinAlpha,
                cosSqAlpha,
                cos2SigmaM,
                C;
        do {
            sinLambda = Math.Sin(lambda);
            cosLambda = Math.Cos(lambda);
            sinSigma = Math.Sqrt((cosU2 * sinLambda) * (cosU2 * sinLambda) + (cosU1 * sinU2 - sinU1 * cosU2 * cosLambda) * (cosU1 * sinU2 - sinU1 * cosU2 * cosLambda));
            if (0 == sinSigma) {
                return 0; // co-incident points
            };
            cosSigma = sinU1 * sinU2 + cosU1 * cosU2 * cosLambda;
            sigma = Math.Atan2(sinSigma, cosSigma);
            sinAlpha = cosU1 * cosU2 * sinLambda / sinSigma;
            cosSqAlpha = 1 - sinAlpha * sinAlpha;
            cos2SigmaM = cosSigma - 2 * sinU1 * sinU2 / cosSqAlpha;
            C = EARTH_BESSEL / 16 * cosSqAlpha * (4 + EARTH_BESSEL * (4 - 3 * cosSqAlpha));
            //  if (isNaN(cos2SigmaM)) {
            //      cos2SigmaM = 0; // equatorial line: cosSqAlpha = 0 (§6)
            //  };
            lambdaP = lambda;
            lambda = L + (1 - C) * EARTH_BESSEL * sinAlpha * (sigma + C * sinSigma * (cos2SigmaM + C * cosSigma * (-1 + 2 * cos2SigmaM * cos2SigmaM)));
        } while (Math.Abs(lambda - lambdaP) > 1e-12 && --iterLimit > 0);

        if (iterLimit == 0) {
            return 0; // formula failed to converge
        };

        double uSq = cosSqAlpha * (EARTH_RADIUS * EARTH_RADIUS - EARTH_RADIUS_MINOR * EARTH_RADIUS_MINOR) / (EARTH_RADIUS_MINOR * EARTH_RADIUS_MINOR),
            A = 1 + uSq / 16384 * (4096 + uSq * (-768 + uSq * (320 - 175 * uSq))),
            B = uSq / 1024 * (256 + uSq * (-128 + uSq * (74 - 47 * uSq))),
            deltaSigma = B * sinSigma * (cos2SigmaM + B / 4 * (cosSigma * (-1 + 2 * cos2SigmaM * cos2SigmaM) - B / 6 * cos2SigmaM * (-3 + 4 * sinSigma * sinSigma) * (-3 + 4 * cos2SigmaM * cos2SigmaM))),
            s = EARTH_RADIUS_MINOR * A * (sigma - deltaSigma);
        return s;
    }

This code was converted from a JavaScript implementation I found here: https://gist.github.com/mathiasbynens/354587

Answer 2

Microsoft may be more correct; but it's not correct!

According the Kallay (2010), SqlGeography STDistance does not return the geodesic distance but the distance along the great ellipse joining the 2 points. The great ellipse distance is not the shortest distance and it does not obey the triangle inequality. The great ellipse returns a reasonable approximation of the geodesic distance for close points. For distant points the error may be as much as 33 km.

An accurate calculation of the geodesic distance is given by the C# library NETGeographicLib (I wrote the underlying C++ library). This is more accurate than Vincenty (10 nanometers instead of 0.5 mm) and, more importantly, the correct result is always obtained. (Vincenty fails to converge for nearly antipodal points.)