2009-08-19 68 views
21

Tôi có một đoạn đường 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.Cách tính 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 cung cấp đ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ả 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.

+2

Bare trong tâm trí rằng Trái Đất, và hệ thống WGS84 thiết kế để xấp xỉ nó, không phải là một quả cầu - vì vậy tính toán dựa trên giả định rằng tôi ở chính xác –

+7

Tôi không biết tại sao mọi người sẽ giả bài tập về nhà .. .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. Trên thực tế, đó là một dự án mini trước đây của tôi ... –

+1

Đôi khi tôi muốn 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

Trả lời

16

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

Tuyên bố từ chối trước tiên. 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 các đầu 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 (có nguồn gốc ở giữa trái đất) - e.g. here.
  2. Tính T, điểm trên đường AB đó là khu vực gần C, sử dụng các sản phẩm 3 vectơ 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 - e.g. here.

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 quan tâm đến khoảng cách giữa C và đoạn đường ngắn hơn, bạn cần phải bước bổ sung để xác minh rằng T thực sự là trên phân đoạn 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.

Tìm điểm gần nhất trên vòng tròn tuyệt vời. Các đầu vào và đầu ra là các mảng dài-2. Các mảng trung gian là chiều 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 đ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 đơn giản kiểm tra nếu điểm T, mà chúng ta biết là trên vòng tròn lớn cùng như A và B, nằm trong phân khúc ngắn hơn của vòng tròn tuyệt vời này. Tuy nhiên có những 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; 
    }  
+1

Có vẻ ổn. Đầ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 độ. –

+0

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

+3

Nó xuất hiện với tôi rằng dòng "trả về (khoảng cách (a, c) aez

2

Hãy thử Distance from a Point to a Great Circle, 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.

0

Khoảng cách ngắn nhất giữa hai điểm trên hình cầu là cạnh nhỏ hơn 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 lập mô hình toán học một cách toán học.

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

+0

nhưng OP không biết điểm thứ 2. Đ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. –

+0

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 đề xuất câu trả lời cho câu hỏi :) – Jimmeh

0

Tôi đang tìm kiếm điều tương tự 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 vòng kết nối tuyệt vời, 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:

This page đề cập đến "Cross-track khoảng cách", mà về cơ bản có vẻ là những gì bạn đang tìm kiếm.

Ngoài ra, trong chuỗi sau 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 khoảng cách đường trên mặt phẳng 2D (với line_locate_point của PostGIS)), và sau đó (2) tính 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 liên kết dưới "có liên quan":

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

1

Đây là mã hoàn chỉnh cho câu trả lời chấp nhận như ideone fiddle (tìm thấy 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(); 

    } 



} 

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 gặp sự cố 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 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 đã đóng.

1

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

Đối với khoảng cách lên đến vài ngàn mét tôi sẽ đơn giản hóa vấn đề từ lĩnh vực để máy bay. Sau đó, vấn đề khá đơn giản vì việc tính toán tam giác dễ dàng có thể được sử dụng:

Chúng tôi có điểm A và B và tìm khoảng cách X đến dò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; 
Các vấn đề liên quan