如何计算球体上从点到线段的距离?

language-agnostic geometry gis geospatial distance

13751 观看

7回复

4889 作者的声誉

我在地球上有一个线段(大圆圈部分)。线段由其末端的坐标定义。显然,两个点定义了两个线段,因此假设我对较短的线段感兴趣。

我给了第三点,我正在寻找线与点之间的(最短)距离。

所有坐标均以经度\纬度(WGS 84)给出。

如何计算距离?

任何合理的编程语言的解决方案都可以。

作者: daphshez 的来源 发布者: 2009 年 8 月 19 日

回应 (7)


3

19368 作者的声誉

尝试数学博士问点到大圆的距离。您仍然需要将经度/纬度转换为球坐标和地球半径的比例,但这似乎是一个很好的方向。

作者: Yuval F 发布者: 19.08.2009 12:43

0

2463 作者的声誉

球面上两个点之间的最短距离是穿过两个点的大圆的较小边。我确定你已经知道了。这里有一个类似的问题http://www.physicsforums.com/archive/index.php/t-178252.html,可以帮助您进行数学建模。

老实说,我不确定您获得此编码示例的可能性有多大。

作者: Jimmeh 发布者: 19.08.2009 12:45

0

12384 作者的声誉

我现在基本上正在寻找相同的东西,只是严格地说我不在乎有一个大圆的段,而是只想到整个圆上任何点的距离。

我目前正在调查两个链接:

本页提到“跨轨距离”,这基本上是您想要的。

另外,在PostGIS邮件列表上的以下线程中,尝试似乎是(1)使用与2D平面上的线距相同的公式(使用PostGIS的line_locate_point)确定大圆上的最近点,然后(2)计算球体上该点与第三个点之间的距离。我不知道数学上的步骤(1)是否正确,但是我会感到惊讶。

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

最后,我只是看到以下链接位于“相关”下:

点到线大圆弧功能的距离不正确。

作者: miracle2k 发布者: 19.08.2009 05:36

18

4889 作者的声誉

决定

这是我自己的解决方案,基于Ask Math博士的想法。很高兴看到您的反馈。

免责声明。该解决方案对于球体是正确的。地球不是球体,坐标系(WGS 84)并不假定它是球体。因此,这只是一个近似值,我无法真正估计出错误。同样,对于很小的距离,也可以通过假设所有事物都只是共面来获得良好的近似值。同样,我不知道距离必须有多小。

现在开始营业。我将线A,B和第三个点C的末端称为。基本上,算法是:

  1. 首先将坐标转换为笛卡尔坐标(原点位于地球中心),例如此处
  2. 使用以下3个矢量积计算T,即AB线上最接近C的点:

    G = A x B

    F = C x G

    T = G x F

  3. 归一化T并乘以地球半径。

  4. 将T转换回经度\纬度。
  5. 计算T和C之间的距离- 例如此处

如果要查找C与A和B定义的大圆之间的距离,这些步骤就足够了。如果像我一样,您对C和较短的线段之间的距离感兴趣,则需要采取额外的步骤来验证T确实在此段上。如果不是,那么最接近的点必然是端点A或B之一-最简单的方法是检查哪一个。

一般而言,以下三个向量乘积的思想如下。第一个(G)给我们A和B大圆的平面(因此,该平面包含A,B和原点)。第二个(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;
} 

这是一种测试点T(我们知道它与A和B在同一大圆上)是否在该大圆的较短段上的简单方法。但是,有更有效的方法可以做到这一点:

   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;
   }    
作者: daphshez 发布者: 19.08.2009 07:53

1

191 作者的声誉

这是可接受的答案的完整代码,如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,这是错误的。

我认为由两点指定的线不是封闭线段。

作者: loleksy 发布者: 11.06.2014 05:27

1

21 作者的声誉

如果有人需要,这是笨拙的答案移植到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;
        }
作者: Dedalo 发布者: 10.10.2015 11:46

1

332 作者的声誉

对于数千米的距离,我将简化从球体到飞机的问题。然后,问题很简单,因为可以使用简单的三角计算:

我们有点A和B,并寻找到线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;
作者: Martin Koubek 发布者: 26.10.2015 10:31
32x32