Câu hỏi Làm thế nào để tính toán khoảng cách từ một điểm đến một đoạn đường thẳng, trên một hình cầu?


Tôi có một đoạn thẳng (phần vòng tròn lớn) trên trái đất. Phân đoạn đường được xác định bởi tọa độ của các đầu của nó. Rõ ràng, hai điểm xác định hai đoạn đường thẳng, vì vậy giả sử tôi quan tâm đến đoạn ngắn hơn.

Tôi được đưa ra một điểm thứ ba, và tôi đang tìm khoảng cách (ngắn nhất) giữa đường thẳng và điểm.

Tất cả các tọa độ được đưa ra theo kinh độ \ vĩ độ (WGS 84).

Làm cách nào để tính khoảng cách?

Một giải pháp trong bất kỳ ngôn ngữ lập trình hợp lý nào sẽ làm.


25
2017-08-19 12:19


gốc


Xin lưu ý rằng hệ thống Trái đất và WGS84 được thiết kế để ước tính nó, không phải là một hình cầu - vì vậy các phép tính dựa trên giả định đó của tôi là không chính xác - John Sibly
Tôi không biết tại sao bất cứ ai sẽ giả định bài tập về nhà ... Tôi đối phó với các điểm về sự gần đúng hình cầu của Trái Đất tại nơi làm việc. Thực tế, đó là một dự án mini trước đây của tôi ... - Thomas Owens
Đôi khi tôi ước tôi vẫn còn ở giai đoạn bài tập về nhà. Điều này tuy nhiên là hoàn toàn làm việc. Nhà vẫn còn khoảng hai giờ nữa. - daphshez
Kinda xúc phạm nhận xét về bài tập về nhà. Người ta sẽ cần ít nhất 2 hoặc 3 học kỳ của chương trình giảng dạy chính toán học từ một trường đại học hàng đầu để xử lý việc này. Có lẽ Lazarus tốt nghiệp MIT từ năm 20 tuổi? - azheglov
Để đơn giản hóa khoảng cách ngắn hơn, hãy xem: stackoverflow.com/questions/20231258/… - Martin Koubek


Các câu trả lời:


Đây là giải pháp của riêng tôi, dựa trên ý tưởng trong hỏi Tiến sĩ Toán. Tôi rất vui khi thấy phản hồi của bạn.

Tuyên bố từ chối trước. Giải pháp này là chính xác cho các lĩnh vực. Trái đất không phải là một hình cầu, và hệ tọa độ (WGS 84) không cho rằng đó là một hình cầu. Vì vậy, đây chỉ là một xấp xỉ, và tôi thực sự không thể ước tính là lỗi. Ngoài ra, với khoảng cách rất nhỏ, có lẽ cũng có thể có được xấp xỉ tốt bằng cách giả định mọi thứ chỉ là một đồng phẳng. Một lần nữa tôi không biết làm thế nào "nhỏ" khoảng cách phải được.

Bây giờ để kinh doanh. Tôi sẽ gọi kết thúc của các dòng A, B và điểm thứ ba C. Về cơ bản, thuật toán là:

  1. chuyển đổi tọa độ đầu tiên thành tọa độ Descartes (với nguồn gốc ở trung tâm của trái đất) - ví dụ. đây.
  2. Tính T, điểm trên đường AB gần nhất với C, sử dụng 3 sản phẩm vector sau:

    G = A x B

    F = C x G

    T = G x F

  3. Chuẩn hóa T và nhân với bán kính của trái đất.

  4. Chuyển đổi T trở lại kinh độ \ vĩ độ.
  5. Tính khoảng cách giữa T và C - ví dụ. đây.

Các bước này là đủ nếu bạn đang tìm khoảng cách giữa C và vòng tròn lớn được xác định bởi A và B. Nếu như tôi bạn quan tâm đến khoảng cách giữa C và đoạn đường ngắn hơn, bạn cần thực hiện thêm bước xác minh T thực sự là trên phân khúc này. Nếu không, thì nhất thiết điểm gần nhất là một trong các đầu A hoặc B - cách dễ nhất là kiểm tra xem cái nào.

Nói chung, ý tưởng đằng sau ba sản phẩm vectơ là như sau. Cái đầu tiên (G) cho chúng ta mặt phẳng của vòng tròn lớn của A và B (vì vậy mặt phẳng chứa A, B và gốc). Thứ hai (F) cho chúng ta vòng tròn lớn đi qua C và vuông góc với G. Sau đó T là giao điểm của các vòng tròn lớn được xác định bởi F và G, được đưa đến độ dài chính xác bằng cách chuẩn hóa và nhân với R.

Đây là một số mã Java một phần để thực hiện nó.

Tìm điểm gần nhất trên vòng tròn lớn. Các đầu vào và đầu ra là các mảng dài-2. Các mảng trung gian có độ dài 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);
}

Tìm điểm gần nhất trên phân đoạn:

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

Đây là một phương pháp thử nghiệm đơn giản nếu điểm T, mà chúng ta biết là trên cùng một vòng tròn lớn như A và B, nằm trên đoạn ngắn hơn của vòng tròn lớn này. Tuy nhiên có nhiều phương pháp hiệu quả hơn để làm điều đó:

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

17
2017-08-19 19:53



Có vẻ tốt. Đầu tiên tôi cũng chuẩn hóa a_, b_ và c_, sao cho mọi thứ ở trên quả bóng đơn vị thay vì trên trái đất. Các sản phẩm vector theo cách này vẫn cung cấp cho bạn các vectơ đơn vị và bạn chỉ phải nhân với bán kính của trái đất để có được các giá trị được chia tỷ lệ chính xác của t và khoảng cách. Ngoài ra, tôi tin rằng nó dễ dàng hơn để tìm khoảng cách trong các tọa độ Descartes (sử dụng định lý Pythagoras) hơn là tìm khoảng cách giữa các vĩ độ kinh độ. - Yuval F
Cảm ơn, chính xác những gì tôi cần. T là trên đoạn mã Minor Arc là những gì tôi đã mất tích. Dave. - daveD
Nó xuất hiện với tôi rằng dòng "trở lại (khoảng cách (a, c) <khoảng cách (b, c))? A: c;" nên được "trả về (khoảng cách (a, c) <khoảng cách (b, c))? a: b;" - aez
tôi ước có ví dụ với số đơn giản - Azam
Một cách tiếp cận khác để xác định xem T là trên đường đi từ A đến B: Giả sử A, B, T là các vector 3d và tất cả chúng đều nằm trong mặt phẳng (Điều này được đảm bảo bởi việc xây dựng T). T là giữa A và B nếu và chỉ khi (AT> AB và BT> AB). Đây * là sản phẩm chấm. Động lực: AT là cosin của góc giữa A và T. Góc giữa A và T cần nhỏ hơn góc giữa A và B. Vì cosin đơn điệu giảm từ 0 xuống Pi, nên dịch sang AT> AB. Quan hệ BT> A * B có động lực tương tự. - thomasfermi


Thử Khoảng cách từ điểm đến vòng kết nối lớn, từ Ask Dr. Math. Bạn vẫn cần phải biến đổi kinh độ / vĩ độ thành các tọa độ hình cầu và tỷ lệ cho bán kính của trái đất, nhưng điều này có vẻ như là một hướng tốt.


3
2017-08-19 12:43





Đây là mã hoàn chỉnh cho câu trả lời được chấp nhận là ideone fiddle (tìm thấy đây):

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();

     }



}

Nó hoạt động tốt cho dữ liệu nhận xét:

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

Nút gần nhất là: 50.17493121381319,19.05846668493702

Nhưng tôi có vấn đề với dữ liệu này:

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

Nút gần nhất là: 52.00834987257176,17.542691313436357 là sai.

Tôi nghĩ rằng dòng được chỉ định bởi hai điểm không phải là một đoạn khép kín.


1
2018-06-11 17:27





Nếu ai đó cần nó đây là câu trả lời loleksy chuyển đến 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;
        }

1
2017-10-10 23:46





Đối với khoảng cách lên đến vài nghìn mét, tôi sẽ đơn giản hóa vấn đề từ hình cầu đến mặt phẳng. Sau đó, vấn đề là khá đơn giản như một phép tính tam giác dễ dàng có thể được sử dụng:

Chúng ta có điểm A và B và tìm kiếm khoảng cách X tới đường AB. Sau đó:

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;

1
2017-10-26 10:31





Khoảng cách ngắn nhất giữa hai điểm trên quả cầu là mặt nhỏ của vòng tròn lớn đi qua hai điểm. Tôi chắc rằng bạn đã biết điều này rồi. Có một câu hỏi tương tự ở đây http://www.physicsforums.com/archive/index.php/t-178252.html có thể giúp bạn mô hình hóa toán học.

Tôi không chắc bạn có khả năng làm thế nào để có được một ví dụ được mã hóa về điều này, thành thật mà nói.


0
2017-08-19 12:45



nhưng OP không biết điểm thứ hai. Điểm P1 = một điểm không nằm trong vòng tròn lớn được xác định bởi đoạn đường, đã biết. Điểm P2 = điểm gần nhất với P1 trên vòng tròn lớn đó, không biết. - Jason S
Tôi hiểu điều đó. Tôi chỉ cần đặt định nghĩa về khoảng cách ngắn nhất giữa hai điểm trên một hình cầu, cùng với một liên kết đến một nơi nào đó hỏi cùng một câu hỏi từ một quan điểm toán học. Tôi không gợi ý câu trả lời cho câu hỏi :) - Jimmeh


Về cơ bản tôi đang tìm kiếm cùng một thứ ngay bây giờ, ngoại trừ việc tôi nghiêm túc nói không quan tâm đến việc có một phân khúc của một vòng tròn lớn, mà chỉ muốn khoảng cách đến bất kỳ điểm nào trên vòng tròn đầy đủ.

Hai liên kết tôi hiện đang điều tra:

Trang này đề cập đến "Khoảng cách theo dõi chéo", về cơ bản có vẻ là những gì bạn đang tìm kiếm.

Ngoài ra, trong chuỗi sau đây trên danh sách gửi thư của PostGIS, nỗ lực dường như (1) xác định điểm gần nhất trên vòng tròn lớn với cùng công thức được sử dụng cho đường thẳng trên mặt phẳng 2D (với line_locate_point của PostGIS), và sau đó (2) tính toán khoảng cách giữa điểm đó và điểm thứ ba trên một hình cầu. Tôi không có ý tưởng nếu bước toán học (1) là chính xác, nhưng tôi sẽ ngạc nhiên.

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

Cuối cùng, tôi chỉ thấy rằng sau đây được liên kết trong "Liên quan":

Khoảng cách từ Point To Line chức năng vòng tròn tuyệt vời không hoạt động đúng.


0