为什么 SqlGeography STDistance 有差异？
Why is there a difference in SqlGeography STDistance?
我一直在追踪我的代码中的一个错误，我发现这是因为 Microsoft c# SqlGeography 2014 库 returns STDistance 的结果与我用于计算之间距离的常规代码略有不同点.
我写了一个小的控制台 exe 来演示这个问题，但我不明白为什么我会得到如此不同的结果？
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 米是 巨大 差异。
好的，经过大量挖掘我找到了答案，Microsoft 更多 正确。
具体来说，他们正在使用 Vincenty's formulae. Accuracy is within 0.5mm (not metre, half a millimetre) instead of 0.3% with Haversine formula。
原因是 Haversine（我用过，Google，显然 Bing 地图太快了）速度很快，但它依赖于球形地球而不是椭球体。 Microsoft 使用椭球体地球而不是球体来计算距离，从而提供更准确的结果。
我像这样在 C# 中实现了 Vincenty 的方法，到目前为止它一直有效，但离生产准备还很远。
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
Microsoft 可能 更多 正确；但不正确！
根据Kallay (2010)，
SqlGeography STDistance 不是 return 测地距离，而是
沿着连接 2 点的大椭圆的距离。最棒的
椭圆距离不是最短距离，不服从
三角不等式大椭圆 return 是合理的
接近点的测地线距离的近似值。对于遥远的
点误差可能高达33公里
测地距离的精确计算由C#给出
图书馆
NETGeographicLib
（我编写了底层 C++ 库）。这比
Vincenty（10 纳米而不是 0.5 毫米），更重要的是，
总能得到正确的结果。 （Vincenty 无法收敛
几乎是对映点。)
我一直在追踪我的代码中的一个错误，我发现这是因为 Microsoft c# SqlGeography 2014 库 returns STDistance 的结果与我用于计算之间距离的常规代码略有不同点.
我写了一个小的控制台 exe 来演示这个问题，但我不明白为什么我会得到如此不同的结果？
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 米是 巨大 差异。
好的，经过大量挖掘我找到了答案，Microsoft 更多 正确。
具体来说，他们正在使用 Vincenty's formulae. Accuracy is within 0.5mm (not metre, half a millimetre) instead of 0.3% with Haversine formula。
原因是 Haversine（我用过，Google，显然 Bing 地图太快了）速度很快，但它依赖于球形地球而不是椭球体。 Microsoft 使用椭球体地球而不是球体来计算距离，从而提供更准确的结果。
我像这样在 C# 中实现了 Vincenty 的方法，到目前为止它一直有效，但离生产准备还很远。
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
Microsoft 可能 更多 正确；但不正确！
根据Kallay (2010)， SqlGeography STDistance 不是 return 测地距离，而是 沿着连接 2 点的大椭圆的距离。最棒的 椭圆距离不是最短距离，不服从 三角不等式大椭圆 return 是合理的 接近点的测地线距离的近似值。对于遥远的 点误差可能高达33公里
测地距离的精确计算由C#给出 图书馆 NETGeographicLib （我编写了底层 C++ 库）。这比 Vincenty（10 纳米而不是 0.5 毫米），更重要的是， 总能得到正确的结果。 （Vincenty 无法收敛 几乎是对映点。)