ポイントからラインセグメントまで、球体上の距離を計算する方法は？

https://stackoverflow.com/questions/1299567

18-09-2019
|

質問

私は地球にラインセグメント（大規模な部分）を持っています。ラインセグメントは、その端の座標によって定義されます。明らかに、2つのポイントが2つのラインセグメントを定義しているため、短いセグメントに興味があると仮定します。

3番目のポイントが与えられており、ラインとポイントの間に（最短）距離を探しています。

すべての座標は、経度緯度（WGS 84）で与えられます。

距離を計算するにはどうすればよいですか？

合理的なプログラミング言語のソリューションが行います。

解決

これが私自身の解決策です。数学博士に聞いてください. 。私はあなたのフィードバックを見てうれしいです。

最初に免責事項。このソリューションは球体に適しています。地球は球体ではなく、座標システム（WGS 84）はそれが球体であるとは想定していません。したがって、これは単なる近似であり、実際にはエラーであるとは推定できません。また、非常にわずかな距離では、すべてが単なる副次的であると仮定することで、良い近似を取得することもおそらく可能です。繰り返しますが、私は距離がどれほど「小さい」かわかりません。

今ビジネスに。行A、B、および3番目のポイントCの端を呼び出します。基本的に、アルゴリズムは次のとおりです。

最初に座標をデカルト座標に変換します（地球の中心に起源があります） - たとえば、ここ.
次の3つのベクトル製品を使用して、Cに最も近いラインABのポイントであるtを計算します。

g = a x b

f = c x g

t = g x f
tを正規化し、地球の半径を掛けます。
Tを経度緯度に戻します。
TとCの間の距離を計算します - たとえば、ここ.

これらの手順は、AとBで定義されたCと大円の間の距離を探している場合に十分です。私のようにCと短いラインセグメントの間の距離に興味がある場合は、それを確認するための追加のステップを実行する必要があります。 Tは確かにこのセグメントにあります。そうでない場合、必ずしも最も近いポイントはAまたはBの端の1つです。最も簡単な方法は、どちらをチェックするかです。

一般的に、3つのベクトル製品の背後にあるアイデアは次のとおりです。最初のもの（g）は、AとBの大規模な円の平面を与えてくれます（したがって、A、B、および原点を含む平面）。 2番目の（f）は、cを通過する大規模な円を与え、Gに垂直です。Tは、fとgによって定義されたグレートサークルの交差点であり、Rによって正規化と乗算によって正しい長さにもたらされます。

これを行うための部分的なJavaコードがいくつかあります。

大規模に最も近いポイントを見つける。入力と出力は長さ2アレイです。中間アレイは長さ3です。

double[] nearestPointGreatCircle(double[] a, double[] b, double c[])
{
    double[] a_ = toCartsian(a);
    double[] b_ = toCartsian(b);
    double[] c_ = toCartsian(c);

    double[] G = vectorProduct(a_, b_);
    double[] F = vectorProduct(c_, G);
    double[] t = vectorProduct(G, F);
    normalize(t);
    multiplyByScalar(t, R_EARTH);
    return fromCartsian(t);
}

セグメントで最も近いポイントを見つける：

double[] nearestPointSegment (double[] a, double[] b, double[] c)
{
   double[] t= nearestPointGreatCircle(a,b,c);
   if (onSegment(a,b,t))
     return t;
   return (distance(a,c) < distance(b,c)) ? a : c;
}

これは、AやBと同じ大規模なサークルにあることがわかっているポイントTが、この大規模なサークルの短いセグメントにある場合、テストする簡単なテスト方法です。ただし、より効率的な方法があります。

   boolean onSegment (double[] a, double[] b, double[] t)
   {
     // should be   return distance(a,t)+distance(b,t)==distance(a,b), 
     // but due to rounding errors, we use: 
     return Math.abs(distance(a,b)-distance(a,t)-distance(b,t)) < PRECISION;
   }

他のヒント

試すポイントから大輪までの距離, 、Ask Dr. Mathから。経度/緯度を球状の座標に変換し、地球の半径のスケールを変える必要がありますが、これは良い方向のようです。

これは、IDEONEフィドルとして受け入れられた回答の完全なコードです（見つかりましたここ):

import java.util.*;
import java.lang.*;
import java.io.*;

/* Name of the class has to be "Main" only if the class is public. */
class Ideone
{



    private static final double _eQuatorialEarthRadius = 6378.1370D;
    private static final double _d2r = (Math.PI / 180D);
    private static double PRECISION = 0.1;





    // Haversine Algorithm
    // source: http://stackoverflow.com/questions/365826/calculate-distance-between-2-gps-coordinates

    private static double HaversineInM(double lat1, double long1, double lat2, double long2) {
        return  (1000D * HaversineInKM(lat1, long1, lat2, long2));
    }

    private static double HaversineInKM(double lat1, double long1, double lat2, double long2) {
        double dlong = (long2 - long1) * _d2r;
        double dlat = (lat2 - lat1) * _d2r;
        double a = Math.pow(Math.sin(dlat / 2D), 2D) + Math.cos(lat1 * _d2r) * Math.cos(lat2 * _d2r)
                * Math.pow(Math.sin(dlong / 2D), 2D);
        double c = 2D * Math.atan2(Math.sqrt(a), Math.sqrt(1D - a));
        double d = _eQuatorialEarthRadius * c;
        return d;
    }

    // Distance between a point and a line

    public static void pointLineDistanceTest() {

        //line
        //double [] a = {50.174315,19.054743};
        //double [] b = {50.176019,19.065042};
        double [] a = {52.00118, 17.53933};
        double [] b = {52.00278, 17.54008};

        //point
        //double [] c = {50.184373,19.054657};
        double [] c = {52.008308, 17.542927};
        double[] nearestNode = nearestPointGreatCircle(a, b, c);
        System.out.println("nearest node: " + Double.toString(nearestNode[0]) + "," + Double.toString(nearestNode[1]));
        double result =  HaversineInM(c[0], c[1], nearestNode[0], nearestNode[1]);
        System.out.println("result: " + Double.toString(result));
    }

    // source: http://stackoverflow.com/questions/1299567/how-to-calculate-distance-from-a-point-to-a-line-segment-on-a-sphere
    private static double[] nearestPointGreatCircle(double[] a, double[] b, double c[])
    {
        double[] a_ = toCartsian(a);
        double[] b_ = toCartsian(b);
        double[] c_ = toCartsian(c);

        double[] G = vectorProduct(a_, b_);
        double[] F = vectorProduct(c_, G);
        double[] t = vectorProduct(G, F);

        return fromCartsian(multiplyByScalar(normalize(t), _eQuatorialEarthRadius));
    }

    @SuppressWarnings("unused")
    private static double[] nearestPointSegment (double[] a, double[] b, double[] c)
    {
       double[] t= nearestPointGreatCircle(a,b,c);
       if (onSegment(a,b,t))
         return t;
       return (HaversineInKM(a[0], a[1], c[0], c[1]) < HaversineInKM(b[0], b[1], c[0], c[1])) ? a : b;
    }

     private static boolean onSegment (double[] a, double[] b, double[] t)
       {
         // should be   return distance(a,t)+distance(b,t)==distance(a,b), 
         // but due to rounding errors, we use: 
         return Math.abs(HaversineInKM(a[0], a[1], b[0], b[1])-HaversineInKM(a[0], a[1], t[0], t[1])-HaversineInKM(b[0], b[1], t[0], t[1])) < PRECISION;
       }


    // source: http://stackoverflow.com/questions/1185408/converting-from-longitude-latitude-to-cartesian-coordinates
    private static double[] toCartsian(double[] coord) {
        double[] result = new double[3];
        result[0] = _eQuatorialEarthRadius * Math.cos(Math.toRadians(coord[0])) * Math.cos(Math.toRadians(coord[1]));
        result[1] = _eQuatorialEarthRadius * Math.cos(Math.toRadians(coord[0])) * Math.sin(Math.toRadians(coord[1]));
        result[2] = _eQuatorialEarthRadius * Math.sin(Math.toRadians(coord[0]));
        return result;
    }

    private static double[] fromCartsian(double[] coord){
        double[] result = new double[2];
        result[0] = Math.toDegrees(Math.asin(coord[2] / _eQuatorialEarthRadius));
        result[1] = Math.toDegrees(Math.atan2(coord[1], coord[0]));

        return result;
    }


    // Basic functions
    private static double[] vectorProduct (double[] a, double[] b){
        double[] result = new double[3];
        result[0] = a[1] * b[2] - a[2] * b[1];
        result[1] = a[2] * b[0] - a[0] * b[2];
        result[2] = a[0] * b[1] - a[1] * b[0];

        return result;
    }

    private static double[] normalize(double[] t) {
        double length = Math.sqrt((t[0] * t[0]) + (t[1] * t[1]) + (t[2] * t[2]));
        double[] result = new double[3];
        result[0] = t[0]/length;
        result[1] = t[1]/length;
        result[2] = t[2]/length;
        return result;
    }

    private static double[] multiplyByScalar(double[] normalize, double k) {
        double[] result = new double[3];
        result[0] = normalize[0]*k;
        result[1] = normalize[1]*k;
        result[2] = normalize[2]*k;
        return result;
    }

     public static void main(String []args){
        System.out.println("Hello World");
        Ideone.pointLineDistanceTest();

     }



}

コメントされたデータには正常に動作します。

//line
double [] a = {50.174315,19.054743};
double [] b = {50.176019,19.065042};
//point
double [] c = {50.184373,19.054657};

最寄りのノードは、50.17493121381319,19.05846668493702です

しかし、私はこのデータに問題があります：

double [] a = {52.00118, 17.53933};
double [] b = {52.00278, 17.54008};
//point
double [] c = {52.008308, 17.542927};

最寄りのノードは、52.00834987257176,17.542691313436357です。

2つのポイントで指定されたラインは閉じたセグメントではないと思います。

誰かがそれを必要とする場合、これはloleksyの回答ですc＃に移植されました

        private static double _eQuatorialEarthRadius = 6378.1370D;
        private static double _d2r = (Math.PI / 180D);
        private static double PRECISION = 0.1;

        // Haversine Algorithm
        // source: http://stackoverflow.com/questions/365826/calculate-distance-between-2-gps-coordinates

        private static double HaversineInM(double lat1, double long1, double lat2, double long2) {
            return  (1000D * HaversineInKM(lat1, long1, lat2, long2));
        }

        private static double HaversineInKM(double lat1, double long1, double lat2, double long2) {
            double dlong = (long2 - long1) * _d2r;
            double dlat = (lat2 - lat1) * _d2r;
            double a = Math.Pow(Math.Sin(dlat / 2D), 2D) + Math.Cos(lat1 * _d2r) * Math.Cos(lat2 * _d2r)
                    * Math.Pow(Math.Sin(dlong / 2D), 2D);
            double c = 2D * Math.Atan2(Math.Sqrt(a), Math.Sqrt(1D - a));
            double d = _eQuatorialEarthRadius * c;
            return d;
        }

        // Distance between a point and a line
        static double pointLineDistanceGEO(double[] a, double[] b, double[] c)
        {

            double[] nearestNode = nearestPointGreatCircle(a, b, c);
            double result = HaversineInKM(c[0], c[1], nearestNode[0], nearestNode[1]);

            return result;
        }

        // source: http://stackoverflow.com/questions/1299567/how-to-calculate-distance-from-a-point-to-a-line-segment-on-a-sphere
        private static double[] nearestPointGreatCircle(double[] a, double[] b, double [] c)
        {
            double[] a_ = toCartsian(a);
            double[] b_ = toCartsian(b);
            double[] c_ = toCartsian(c);

            double[] G = vectorProduct(a_, b_);
            double[] F = vectorProduct(c_, G);
            double[] t = vectorProduct(G, F);

            return fromCartsian(multiplyByScalar(normalize(t), _eQuatorialEarthRadius));
        }

        private static double[] nearestPointSegment (double[] a, double[] b, double[] c)
        {
           double[] t= nearestPointGreatCircle(a,b,c);
           if (onSegment(a,b,t))
             return t;
           return (HaversineInKM(a[0], a[1], c[0], c[1]) < HaversineInKM(b[0], b[1], c[0], c[1])) ? a : b;
        }

         private static bool onSegment (double[] a, double[] b, double[] t)
           {
             // should be   return distance(a,t)+distance(b,t)==distance(a,b), 
             // but due to rounding errors, we use: 
             return Math.Abs(HaversineInKM(a[0], a[1], b[0], b[1])-HaversineInKM(a[0], a[1], t[0], t[1])-HaversineInKM(b[0], b[1], t[0], t[1])) < PRECISION;
           }


        // source: http://stackoverflow.com/questions/1185408/converting-from-longitude-latitude-to-cartesian-coordinates
        private static double[] toCartsian(double[] coord) {
            double[] result = new double[3];
            result[0] = _eQuatorialEarthRadius * Math.Cos(deg2rad(coord[0])) * Math.Cos(deg2rad(coord[1]));
            result[1] = _eQuatorialEarthRadius * Math.Cos(deg2rad(coord[0])) * Math.Sin(deg2rad(coord[1]));
            result[2] = _eQuatorialEarthRadius * Math.Sin(deg2rad(coord[0]));
            return result;
        }

        private static double[] fromCartsian(double[] coord){
            double[] result = new double[2];
            result[0] = rad2deg(Math.Asin(coord[2] / _eQuatorialEarthRadius));
            result[1] = rad2deg(Math.Atan2(coord[1], coord[0]));

            return result;
        }


        // Basic functions
        private static double[] vectorProduct (double[] a, double[] b){
            double[] result = new double[3];
            result[0] = a[1] * b[2] - a[2] * b[1];
            result[1] = a[2] * b[0] - a[0] * b[2];
            result[2] = a[0] * b[1] - a[1] * b[0];

            return result;
        }

        private static double[] normalize(double[] t) {
            double length = Math.Sqrt((t[0] * t[0]) + (t[1] * t[1]) + (t[2] * t[2]));
            double[] result = new double[3];
            result[0] = t[0]/length;
            result[1] = t[1]/length;
            result[2] = t[2]/length;
            return result;
        }

        private static double[] multiplyByScalar(double[] normalize, double k) {
            double[] result = new double[3];
            result[0] = normalize[0]*k;
            result[1] = normalize[1]*k;
            result[2] = normalize[2]*k;
            return result;
        }

数千メートルまでの距離については、球体から飛行機までの問題を簡素化します。次に、簡単な三角形の計算を使用できるので、問題はかなり単純です。

ポイントAとBがあり、ABからABへの距離xを探します。それで：

Location a;
Location b;
Location x;

double ax = a.distanceTo(x);
double alfa = (Math.abs(a.bearingTo(b) - a.bearingTo(x))) / 180
            * Math.PI;
double distance = Math.sin(alfa) * ax;

球上の2つのポイント間の最短距離は、2つのポイントを通る大円の小さな側面です。私はあなたがすでにこれを知っていると確信しています。ここにも同様の質問があります http://www.physicsforums.com/archive/index.php/t-178252.html それはあなたがそれを数学的にモデル化するのに役立つかもしれません。

正直に言うと、このコード化された例を取得する可能性がどれほどかかるかはわかりません。

私は基本的に同じことを探していますが、私は厳密に言えば、大規模な円のセグメントを持つことを気にしないでください。

私が現在調査している2つのリンク：

このページ「クロストラック距離」に言及していますが、これは基本的にあなたが探しているもののようです。

また、PostGISメーリングリストの次のスレッドでは、（1）2Dプレーン（PostGISのline_locate_pointを使用）のライン距離に使用されるのと同じ式で、大規模に最も近いポイントを決定し、次に試みます。（2）スフェロイド上の3ポイントと3番目のポイントの間の距離を計算します。数学的にステップ（1）が正しいかどうかはわかりませんが、驚かれることになります。

http://postgis.refractions.net/pipermail/postgis-users/2009-july/023903.html

最後に、「関連」の下に次のものがリンクされていることがわかりました。

ポイントからラインへの距離グレートサークル機能は正しく機能しません。

ライセンス： CC-BY-SA と帰属

所属していません StackOverflow