C#：経度緯度の２点間の距離を求めるには

以前、クラスライブラリを使用して距離を取得しましたが、そのときのライブラリが今では使いづらそうなので、算出するコードを作成してみました。
ChatGPT で生成して、少し直した程度です。凄い時代になりましたね。
Vincenty 法で求めると精度が高いようなので、それで生成させました。

using System;

var location1 = new Location(35.6895, 139.6917); // 東京
var location2 = new Location(34.6937, 135.5022); // 大阪

var distance = Vincenty.Distance(location1, location2);
Console.WriteLine($"Distance between Tokyo and Osaka: {distance} m");

public class Location
{
    public double Latitude { get; set; }
    public double Longitude { get; set; }

    public Location(double latitude, double longitude)
    {
        Latitude = latitude;
        Longitude = longitude;
    }
}

public static class Vincenty
{
    private const double Flattening = 1 / 298.257223563; // WGS84楕円体の扁平率
    private const double EquatorialRadius = 6378137.0; // WGS84楕円体の赤道半径

    public static double Distance(Location loc1, Location loc2)
    {
        var f = Flattening;
        var a = EquatorialRadius;
        var b = a * (1 - f);
        var latRad1 = loc1.Latitude * Math.PI / 180.0;
        var lonRad1 = loc1.Longitude * Math.PI / 180.0;
        var latRad2 = loc2.Latitude * Math.PI / 180.0;
        var lonRad2 = loc2.Longitude * Math.PI / 180.0;

        var L = lonRad2 - lonRad1;
        var tanU1 = (1 - f) * Math.Tan(latRad1);
        var cosU1 = 1 / Math.Sqrt((1 + tanU1 * tanU1));
        var sinU1 = tanU1 * cosU1;
        var tanU2 = (1 - f) * Math.Tan(latRad2);
        var cosU2 = 1 / Math.Sqrt((1 + tanU2 * tanU2));
        var sinU2 = tanU2 * cosU2;

        double lambda = L;
        double lambdaP = 2 * Math.PI;

        double sinLambda = 0;
        double cosLambda = 0;
        double sinSigma = 0;
        double cosSigma = 0;
        double sigma = 0;
        double sinAlpha = 0;
        double cosSqAlpha = 0;
        double cos2SigmaM = 0;
        double C = 0;

        while (Math.Abs(lambda - lambdaP) > 1e-12)
        {
            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 (sinSigma == 0)
            {
                return 0;  // 2点間の距離が0の場合
            }

            cosSigma = sinU1 * sinU2 + cosU1 * cosU2 * cosLambda;
            sigma = Math.Atan2(sinSigma, cosSigma);
            sinAlpha = cosU1 * cosU2 * sinLambda / sinSigma;
            cosSqAlpha = 1 - sinAlpha * sinAlpha;
            if (cosSqAlpha == 0)
            {
                cos2SigmaM = 0;
            }
            else
            {
                cos2SigmaM = cosSigma - 2 * sinU1 * sinU2 / cosSqAlpha;
            }

            C = f / 16 * cosSqAlpha * (4 + f * (4 - 3 * cosSqAlpha));
            lambdaP = lambda;
            lambda = L + (1 - C) * f * sinAlpha * (sigma + C * sinSigma *
                      (Math.Cos(2 * sigma) + C * cosSigma * (-1 + 2 * Math.Cos(2 * sigma))));
        }
        var uSq = cosSqAlpha * (a * a - b * b) / (b * b);
        var A = 1 + uSq / 16384 * (4096 + uSq * (-768 + uSq * (320 - 175 * uSq)));
        var B = uSq / 1024 * (256 + uSq * (-128 + uSq * (74 - 47 * uSq)));
        var deltaSigma = B * sinSigma * (cos2SigmaM + B / 4 * (cosSigma * (-1 + 2 * cos2SigmaM * cos2SigmaM) -
                         B / 6 * cos2SigmaM * (-3 + 4 * sinSigma * sinSigma) * (-3 + 4 * cos2SigmaM * cos2SigmaM)));
        var s = b * A * (sigma - deltaSigma);

        return s;
    }
}

Distance between Tokyo and Osaka : 397192.2878195703 m

コメント