Wie berechnet man den Abstand von einem Punkt zu einem Liniensegment auf einer Kugel?
https://stackoverflow.com/questions/1299567
-
18-09-2019 - |
italiano
english
français
española
中国
日本の
العربية
Deutsch
한국어
Português
RussianFrage
Ich habe ein Liniensegment (Großkreis -Teil) auf der Erde. Das Liniensegment wird durch die Koordinaten seiner Enden definiert. Offensichtlich definieren zwei Punkte zwei Liniensegmente. Nehmen wir also an, ich interessiere mich für die kürzere.
Ich habe einen dritten Punkt und suche nach dem (kürzesten) Abstand zwischen der Linie und dem Punkt.
Alle Koordinaten sind in Längengrad Breitengrad (WGS 84) angegeben.
Wie berechnet ich die Entfernung?
Eine Lösung in einer angemessenen Programmiersprache ist.
Lösung
Hier ist meine eigene Lösung, basierend auf der Idee in Fragen Sie Dr. Math. Ich würde mich freuen, Ihr Feedback zu sehen.
Haftungsausschluss zuerst. Diese Lösung ist für Kugeln korrekt. Die Erde ist keine Sphäre, und das Koordinatensystem (WGS 84) geht nicht davon aus, dass es sich um eine Kugel handelt. Das ist also nur eine Annäherung, und ich kann nicht wirklich schätzen, dass er ein Fehler ist. Für sehr kleine Entfernungen ist es wahrscheinlich auch möglich, eine gute Annäherung zu erhalten, indem angenommen wird, dass alles nur ein Coplanar ist. Ich weiß wieder nicht, wie "klein" die Entfernungen sein müssen.
Jetzt zum Geschäft. Ich werde die Enden der Zeilen A, B und den dritten Punkt C nennen. Grundsätzlich ist der Algorithmus:
- Konvertieren Sie die Koordinaten zuerst in kartesische Koordinaten (mit dem Ursprung im Mittelpunkt) - zB hier.
Berechnen Sie t den Punkt auf der Linie AB, der C am nächsten liegt, unter Verwendung der folgenden 3 Vektorprodukte:
G = a x b
F = c x g
T = g x f
Normalisieren Sie T und multiplizieren Sie sie mit dem Radius der Erde.
- Konvertieren Sie T wieder in den Breitengrad.
- Berechnen Sie den Abstand zwischen t und c - zB hier.
Diese Schritte reichen aus, wenn Sie nach dem Abstand zwischen C und dem von A und B definierten großen Kreis suchen T ist in der Tat in diesem Segment. Wenn dies nicht der Fall ist, ist der nächste Punkt notwendigerweise einer der Enden A oder B - der einfachste Weg, um zu überprüfen, welche.
Im Allgemeinen ist die Idee hinter den drei Vektorprodukten die folgende. Der erste (g) gibt uns die Ebene des großen Kreises von A und B (so die Ebene, die A, B und den Ursprung enthält). Der zweite (f) gibt uns den großen Kreis, der durch C geht und senkrecht zum G. Dann ist t der Schnittpunkt der großen Kreise, die durch F und G definiert werden, die durch Normalisierung und Multiplikation durch R auf die richtige Länge gebracht werden.
Hier ist ein teilweise Java -Code dafür.
Finden Sie den nächsten Punkt im Großen Kreis. Die Eingänge und Ausgabe sind Länge-2-Arrays. Die Zwischenarrays sind von Länge 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);
}
Finden Sie den nächsten Punkt auf dem Segment:
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;
}
Dies ist eine einfache Testmethode zum Testen, wenn Punkt T, von dem wir wissen, dass es sich um den gleichen Großkreis wie A und B befindet, auf dem kürzeren Segment dieses großen Kreises. Es gibt jedoch effizientere Methoden, um dies zu tun:
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;
}
Andere Tipps
Versuchen Entfernung von einem Punkt zu einem großen Kreis, von Ask Dr. Math. Sie müssen immer noch den Längengrad/den Breitengrad in kugelförmige Koordinaten und die Skala für den Erdradius verwandeln, aber dies scheint eine gute Richtung zu sein.
Dies ist ein vollständiger Code für die akzeptierte Antwort als ideone Fiddle (gefunden hier):
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();
}
}
Es funktioniert gut für kommentierte Daten:
//line
double [] a = {50.174315,19.054743};
double [] b = {50.176019,19.065042};
//point
double [] c = {50.184373,19.054657};
Der nächste Knoten ist: 50.17493121381319, 19905846668493702
Aber ich habe ein Problem mit diesen Daten:
double [] a = {52.00118, 17.53933};
double [] b = {52.00278, 17.54008};
//point
double [] c = {52.008308, 17.542927};
Der nächste Knoten ist: 52.00834987257176,17.542691313436357, was falsch ist.
Ich denke, diese durch zwei Punkte angegebene Linie ist kein geschlossenes Segment.
Wenn jemand es braucht, ist dies die Loleksy -Antwort auf C# portiert
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;
}
Für die Entfernung bis zu ein paar tausend Metern würde ich das Problem von der Kugel zu Ebene vereinfachen. Dann ist das Problem ziemlich einfach, da eine einfache Dreiecksberechnung verwendet werden kann:
Wir haben Punkte A und B und suchen nach einer Entfernung X, um AB zu leiten. Dann:
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;
Der kürzeste Abstand zwischen zwei Punkten auf einer Kugel ist die kleinere Seite des großen Kreises durch die beiden Punkte. Ich bin sicher, Sie wissen das schon. Hier gibt es eine ähnliche Frage http://www.physicsforums.com/archive/index.php/t-178252.html Das kann Ihnen helfen, es mathmatisch zu modellieren.
Ich bin mir nicht sicher, wie wahrscheinlich es ist, dass Sie ein codiertes Beispiel dafür erhalten, um ehrlich zu sein.
Im Grunde suche ich im Grunde nach demselben Ding, außer dass es mir streng genommen egal ist, ein Segment eines großen Kreises zu haben, sondern nur die Entfernung zu einem beliebigen Punkt auf dem Kreis.
Zwei Links, die ich derzeit untersuche:
Diese Seite Erwähnt "Cross-Track Distance", was im Grunde genommen das sein scheint, was Sie suchen.
Außerdem scheint der Versuch im folgenden Thread auf der Postgis-Mailingliste (1) den engsten Punkt auf dem großen Kreis mit der gleichen Formel zu bestimmen, die für die Zeilendistanz auf einer 2D-Ebene verwendet wird (mit postgis line_locate_point), und dann (2) Berechnung des Abstands zwischen diesem und dem dritten Punkt auf einem Sphäroid. Ich habe keine Ahnung, ob mathematisch Schritt (1) korrekt ist, aber ich wäre überrascht.
http://postgis.refractions.net/pipermail/postgis-users/2009-july/023903.html
Schließlich habe ich gerade gesehen, dass das folgende unter "verwandte" verknüpft ist:
Entfernung von Punkt zur Linie Großkreis funktioniert nicht richtig.
