为什么SqlGeography STDistance有区别? [英] Why is there a difference in SqlGeography STDistance?
问题描述
我一直在跟踪代码中的错误，我发现这是因为Microsoft c#SqlGeography 2014库返回的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 regular code for calculating the distance between points.
我编写了一个小型控制台exe来演示该问题，但我不知道为什么会得到如此不同的结果?
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
}
Microsoft是否使用其他计算方法?计算不到1.5公里的距离时，距离超过1米是一个巨大差异.
Is Microsoft using a different calculation method? Being off by over 1 meter when calculating distances less than 1.5Km is a huge discrepancy.
推荐答案
好的，所以经过大量挖掘，我找到了答案，微软更多是正确的.
Ok, so after much digging I found the answer, and Microsoft is more correct.
具体来说，他们使用的是 Vincenty的公式.精度在0.5毫米以内(不是米，半毫米)，而不是
Specifically, they are using Vincenty's formulae. Accuracy is within 0.5mm (not metre, half a millimetre) instead of 0.3% with Haversine formula.
原因是Haversine(我，Google和Bing Maps显然使用过)速度很快，但依赖于球形地球而不是椭圆形. Microsoft使用椭球体Earth而不是球体来计算距离，从而提供了更准确的结果.
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.
我像这样在c#中实现了Vincenty的方法，到目前为止，它已经可以使用了，但是还没准备好投入生产.
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;
}
此代码是从我在此处找到的JavaScript实现转换而成的: https://gist.github.com/mathiasbynens/354587
This code was converted from a JavaScript implementation I found here: https://gist.github.com/mathiasbynens/354587
这篇关于为什么SqlGeography STDistance有区别?的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
