Why Is There A Difference In SqlGeography STDistance?

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 r

Solution 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

Solution 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.)