2009-08-19 16 views
21

Mam odcinek linii (świetny okrąg) na ziemi. Segment linii jest zdefiniowany przez współrzędne jego końców. Oczywiście dwa punkty definiują dwa odcinki linii, więc przyjmijmy, że jestem zainteresowany tym krótszym.Jak obliczyć odległość od punktu do segmentu linii, na kuli?

Dostałem trzeci punkt i szukam (najkrótszej) odległości między linią a punktem.

Wszystkie współrzędne podane są w długości geograficznej \ szerokości geograficznej (WGS 84).

Jak obliczyć odległość?

Rozwiązanie w dowolnym rozsądnym języku programowania.

+2

Bare na uwadze, że i układ WGS84 Ziemia zaprojektowany, aby go zbliżyć, nie jest kulą - tak Obliczenia oparte na założeniu, że moje niedokładne –

+7

nie wiem, dlaczego ktoś miałby zakładać domową .. Zajmuję się punktami na sferycznej aproksymacji Ziemi w pracy. W rzeczywistości był to mój poprzedni mini-projekt ... –

+1

Czasami żałuję, że nie byłem jeszcze na etapie prac domowych. To jednak jest całkowicie praca. Dom jest nadal oddalony o dwie godziny. – daphshez

Odpowiedz

16

Oto moje własne rozwiązanie, oparte na idei w ask Dr. Math. Byłbym szczęśliwy widząc twoją opinię.

Zastrzeżenie pierwsze. To rozwiązanie jest poprawne w przypadku sfer. Ziemia nie jest kulą, a układ współrzędnych (WGS 84) nie zakłada, że ​​jest to sfera. Jest to tylko przybliżenie i nie mogę naprawdę oszacować błędu. Ponadto, w przypadku bardzo małych odległości, możliwe jest również uzyskanie dobrego przybliżenia, zakładając, że wszystko jest po prostu współpłaszczyznowe. Znowu nie wiem, jak "małe" muszą być odległości.

Teraz do firmy. Nazwiemy końce linii A, B i trzeci punkt C.Zasadniczo, algorytm jest:

  1. konwersji współrzędnych kartezjańskich współrzędnych pierwszy (z początku w środku ziemi) - e.g. here.
  2. Oblicz T punkt na linii AB, który jest najbliższy ° C, stosując następujące produkty: 3 wektor:

    G = A x B

    C = C x g

    T = G x F

  3. Normalizuj T i pomnóż przez promień ziemi.

  4. Konwertuj T z powrotem na długość geograficzną \ szerokość geograficzną.
  5. Oblicz odległość między T a C - e.g. here.

Te kroki wystarczą, jeśli szukasz odległości między C a wielkim okręgiem zdefiniowanym przez A i B. Jeśli tak jak ja interesuje cię odległość między C i krótszym odcinkiem linii, musisz wziąć dodatkowy krok polegający na sprawdzeniu, czy T rzeczywiście należy do tego segmentu. Jeśli nie jest, to koniecznie najbliższym punktem jest jeden z końców A lub B - najłatwiejszym sposobem jest sprawdzenie, który z nich.

W kategoriach ogólnych, idea trzech produktów wektorowych jest następująca. Pierwszy (G) daje nam płaszczyznę wielkiego koła A i B (tak więc płaszczyzna zawierająca A, B i pochodzenie). Drugi (F) daje nam wielkie koło, które przechodzi przez C i jest prostopadłe do G. Następnie T jest przecięciem wielkich okręgów zdefiniowanych przez F i G, doprowadzonych do właściwej długości przez normalizację i pomnożenie przez R.

Oto niektóre częściowy kod Java do robienia tego.

Znajdowanie najbliższego punktu w wielkim kręgu. Wejścia i wyjścia są tablicami długości 2. Macierze pośrednie są od długości 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); 
} 

Znalezienie najbliższego punktu na segmencie:

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

Jest to prosta metoda badawcza, jeżeli punkt T, co wiemy, jest na tym samym wielkim kole jako A i B znajduje się na krótszym odcinku tego wielkiego koła. Jednakże istnieją bardziej skuteczne sposoby aby to zrobić:

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

Wygląda dobrze. Najpierw też normalizowałbym a_, b_ i c_, aby wszystko było na kuli jednostki zamiast na ziemi. W ten sposób produkty wektorowe nadal dostarczają wektorów jednostkowych, a ty musisz tylko pomnożyć przez promień Ziemi, aby uzyskać poprawnie skalowane wartości t oraz odległość. Uważam też, że łatwiej jest znaleźć odległości we współrzędnych kartezjańskich (używając twierdzenia Pitagorasa) niż znaleźć odległości między punktami szerokości i długości geograficznej. –

+0

Dzięki, dokładnie to, czego potrzebowałem. T było na mniejszym odcinku kodu łuku było tym, czego mi brakowało. Dave. – daveD

+3

Wydaje mi się, że linia "return (odległość (a, c) aez

2

Wypróbuj Distance from a Point to a Great Circle, od Ask Dr. Math. Nadal musisz przekształcić długość/szerokość geograficzną na współrzędne sferyczne i skalę dla promienia Ziemi, ale wydaje się to dobrym kierunkiem.

0

Najkrótsza odległość między dwoma punktami na kuli to mniejsza strona wielkiego koła przechodzącego przez dwa punkty. Jestem pewien, że już to wiesz. Istnieje podobne pytanie tutaj: http://www.physicsforums.com/archive/index.php/t-178252.html, które może pomóc w modelowaniu go matematycznie.

Nie jestem pewien, jak bardzo prawdopodobne jest, że otrzymasz zakodowany przykład tego, szczerze mówiąc.

+0

, ale OP nie znać drugi punkt. Punkt P1 = punkt nie na wielkim okręgu określonym przez odcinek linii, znany. Punkt P2 = najbliższy punkt P1 na tym wielkim okręgu, nieznany. –

+0

Rozumiem to. Po prostu umieszczam definicję najkrótszej odległości między dwoma punktami na kuli, wraz z linkiem do czegoś zadającego to samo pytanie z matematycznego punktu widzenia. Nie sugeruję odpowiedzi na pytanie :) – Jimmeh

0

Zasadniczo szukam teraz tego samego, z tym wyjątkiem, że ściśle mówiąc nie dbam o to, aby mieć segment wielkiego koła, ale po prostu chcę odległości do dowolnego punktu w pełnym kole.

dwa linki Jestem obecnie dochodzenie:

This page wspomina „Cross-track odległość”, która w zasadzie wydaje się być to, czego szukasz.

Ponadto, w następującym wątku na liście mailingowej PostGIS, próba (1) określa najbliższy punkt na wielkim okręgu z tą samą formułą, używaną dla odległości linii na płaszczyźnie 2D (z PostGIS 'line_locate_point), a następnie (2) obliczenie odległości między tym a trzecim punktem na sferoidzie. Nie mam pojęcia, czy matematycznie krok (1) jest poprawny, ale byłbym zaskoczony.

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

Wreszcie, po prostu zobaczył, że następuje powiązany w sekcji "Powiązane":

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

1

Jest to kompletny kod dla akceptowanych odpowiedzi jako ideone skrzypce (znaleziono 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(); 

    } 



} 

to działa dobrze dla skomentował dane:

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

Najbliższy węzeł: 50.17493121381319,19.05846668493702

Ale mam problem z tymi danymi:

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

Najbliższy węzeł jest: 52.00834987257176,17.542691313436357 co jest źle.

Myślę, że linia określona przez dwa punkty nie jest zamkniętym segmentem.

1

Jeśli ktoś potrzebuje to odpowiedź loleksy przeniesiony do 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

na odległość do kilku tysięcy metrów chciałbym uprościć problem z kuli na płaszczyźnie. Następnie problem jest prosty, ponieważ można użyć prostego obliczenia trójkąta:

Mamy punkty A i B i szukamy odległości X do linii AB. Następnie:

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; 
Powiązane problemy