【问题标题】:How to calculate distance from a point to a line segment, on a sphere?如何计算球体上点到线段的距离?
【发布时间】:2010-11-20 22:30:00
【问题描述】:

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

我得到了第三个点,我正在寻找线和点之间的(最短)距离。

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

如何计算距离?

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

【问题讨论】:

  • 请记住,地球和设计用于近似它的 WGS84 系统不是球体 - 因此基于该假设的计算不准确
  • 我不知道为什么有人会承担家庭作业...我在工作中处理地球近似球形上的点。其实这是我之前的一个小项目……
  • 我有时真希望我还在做作业。然而,这完全是工作。离家还有大约两个小时的路程。
  • 对家庭作业的侮辱性评论。一个人需要至少 2 或 3 个学期的一流大学的数学专业课程来处理这个问题。也许拉撒路在 20 岁时从麻省理工学院毕业?
  • 为了更短的距离做简化,见:stackoverflow.com/questions/20231258/…

标签: language-agnostic geometry gis geospatial distance


【解决方案1】:

这是我自己的解决方案,基于ask Dr. Math 中的想法。我很高兴看到您的反馈。

首先声明。此解决方案适用于球体。地球不是球体,坐标系 (WGS 84) 也不假定它是球体。所以这只是一个近似值,我无法真正估计是错误。此外,对于非常小的距离,通过假设一切都是共面的,也可能获得很好的近似值。同样,我不知道距离必须有多“小”。

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

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

    G = A x B

    F = C x G

    T = G x F

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

  4. 将 T 转换回经度\纬度。
  5. 计算 T 和 C 之间的距离 - e.g. here

如果您正在寻找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;
   }    

【讨论】:

  • 看起来不错。我还将首先对 a_、b_ 和 c_ 进行归一化,以便所有内容都在单位球上而不是地球上。这种方式矢量产品仍然为您提供单位矢量,您只需乘以地球的半径即可获得正确缩放的 t 值和距离。另外,我相信在笛卡尔坐标中找到距离(使用毕达哥拉斯定理)比找到经纬度点之间的距离更容易。
  • 谢谢,正是我需要的。 T is on Minor Arc 段代码是我所缺少的。戴夫。
  • 在我看来,“return (distance(a,c)
  • 我希望有简单数字的例子
  • 另一种判断 T 是否在从 A 到 B 的路径的方法:假设 A,B,T 是 3d 向量并且它们都在一个平面上(这由 T 的构造来保证) . T 在 A 和 B 之间当且仅当 (AT > AB 和 BT > AB)。这里 * 是点积。动机:AT 是 A 和 T 之间角度的余弦。A 和 T 之间的角度需要小于 A 和 B 之间的角度。由于余弦从 0 到 Pi 单调递减,因此转换为 A T > AB。关系 BT > A*B 的动机类似。
【解决方案2】:

尝试来自 Ask Dr. Math 的 Distance from a Point to a Great Circle。您仍然需要将经度/纬度转换为球坐标和地球半径的比例,但这似乎是一个不错的方向。

【讨论】:

    【解决方案3】:

    如果有人需要,这是移植到 c# 的 loleksy 答案

            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;
            }
    

    【讨论】:

      【解决方案4】:

      我现在基本上在寻找相同的东西,除了严格来说我不关心有一个大圆的一部分,而只是想要到整个圆上任何一点的距离。

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

      This page 提到“Cross-track distance”,这基本上就是你要找的。​​p>

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

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

      最后,我刚刚在“相关”下看到以下链接:

      Distance from Point To Line great circle function not working right.

      【讨论】:

        【解决方案5】:

        这是作为ideone fiddle接受答案的完整代码(找到here):

        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 这是错误的。

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

        【讨论】:

          【解决方案6】:

          对于长达几千米的距离,我会简化从球体到平面的问题。 然后,问题就很简单了,因为可以使用简单的三角形计算:

          我们有点 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;
          

          【讨论】:

            【解决方案7】:

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

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

            【讨论】:

            • 但 OP 不知道第二点。点 P1 = 不在由线段定义的大圆上的点,已知。点 P2 = 在那个大圆上离 P1 最近的点,未知。
            • 我明白这一点。我只是简单地定义了球体上两点之间的最短距离,以及从数学角度提出相同问题的某个地方的链接。我不是建议回答这个问题:)
            • 我知道这是 11 年前的事了,所以,@Jimmeh,请忽略这条评论。但对于未来阅读本文的人:这类事情应该是评论,而不是答案,因为它不会尝试回答问题。
            猜你喜欢
            • 1970-01-01
            • 1970-01-01
            • 1970-01-01
            • 1970-01-01
            • 2013-01-26
            • 2018-02-15
            • 1970-01-01
            • 1970-01-01
            • 2020-12-08
            相关资源
            最近更新 更多