2017-09-11 6 views
1

私は自分のコードのバグを追跡していますが、Microsoftのc#SqlGeography 2014ライブラリが、STDistanceのポイントとポイントの距離を計算するための私の通常のコードとは少し異なる結果を返すからです。なぜSqlGeography 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メートル以上離れていると巨大なの矛盾があります。

答えて

2

よく分かりましたので、私は答えを見つけたので、マイクロソフトはより多くです。

具体的には、Vincenty's formulaeを使用しています。精度はHaversine formulaで0.3%の代わりに0.5mm以内です(メートル,半ミリメートル)。

なぜなら、Haversine(私、Google、Bing Mapsがあまりにも明らかに使用した)は高速ですが、楕円の代わりに球状の地球に依存しているからです。マイクロソフトでは、より正確な結果を提供する球の代わりに距離を計算するために楕円地球を使用しています。

私はこのように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

0

Microsoftは正しい以上であってもよいです。それはではない正しい!

Kallay (2010)によれば、 SqlGeography STDistanceは測地線距離が、2点を結ぶ大きな楕円に沿っ 距離を返しません。 楕円距離は最短距離ではなく、 三角不等式には従いません。大きな楕円は、近い点の測地線距離の近似値である を返します。遠方の ポイントの場合、エラーは33kmもある可能性があります。

測地距離の正確な計算

は、C# ライブラリ NETGeographicLib で与えられる(Iは、基礎となるC++ライブラリを書いています)。これは、 Vincenty(0.5ナノメートルの代わりに10ナノメートル)よりも正確であり、さらに重要なことに、正確な結果は常に得られる。 (Vincentyは にほぼ収束しません)

関連する問題