Cómo calcular la distancia de un punto a un segmento de línea, en una esfera?
https://stackoverflow.com/questions/1299567
-
18-09-2019 - |
italiano
english
français
española
中国
日本の
العربية
Deutsch
한국어
Português
RussianPregunta
Tengo un segmento de línea (gran parte círculo) en la tierra. El segmento de línea se define por las coordenadas de sus extremos. Obviamente, dos puntos definen dos segmentos de línea, por lo que asumen Estoy interesado en el más corto.
Me dan un tercer punto, y estoy buscando para la distancia (corta) entre la línea y el punto.
Todas las coordenadas se dan en la longitud \ latitud (WGS 84).
¿Cómo se calcula la distancia?
Una solución en cualquier lenguaje de programación razonable va a hacer.
Solución
Esta es mi propia solución, basada en la idea de rel="nofollow preguntarle al Dr. Matemáticas . Yo estaría feliz de ver sus comentarios.
Negación primera. Esta solución es correcta para esferas. La Tierra no es una esfera, y el sistema de coordenadas (WGS 84) no asuma que es una esfera. Así que esto es sólo una aproximación, y no puedo estimar que es error. Además, para distancias muy pequeñas, es probable que también sea posible para obtener una buena aproximación al suponer que todo es un solo un mismo plano. Una vez más no sé qué "pequeña" las distancias tienen que ser.
Ahora al grano. Voy a llamar a los extremos de las líneas A, B y el tercer punto C. Básicamente, el algoritmo es:
- convertir las coordenadas primero en coordenadas cartesianas (con el origen en el centro de la tierra) - ejemplo, aquí.
-
Calcular T, el punto de la línea AB que es la más cercana a C, usando los productos vectoriales siguiente 3:
G = A x B
F = C x G
T = G x F
-
Normalizar T y se multiplica por el radio de la tierra.
- Convertir T de nuevo a la longitud \ latitud.
- Calcular la distancia entre T y C - por ejemplo aquí.
Estos pasos son suficientes si se busca la distancia entre C y el gran círculo definido por A y B. Si como yo, está interesado en la distancia entre C y el segmento de línea más corto, es necesario tomar la medida adicional de verificar que T es de hecho en este segmento. Si no es así, entonces necesariamente el punto más cercano es uno de los extremos A o B - la forma más fácil es comprobar cuál.
En términos generales, la idea detrás de los tres productos vectoriales es el siguiente. La primera (G) nos da la el plano de la gran círculo de A y B (por lo que el plano que contiene A, B y el origen). La segunda (F) nos da el gran círculo la pasa a través de C y es perpendicular a la G. Entonces T es la intersección de los grandes círculos definidos por F y G, traído a la longitud correcta por la normalización y la multiplicación por R.
Aquí hay algo de código Java parcial para hacerlo.
Encontrar el punto en el círculo más cercano. Las entradas y salida son de longitud-2 arrays. Los arrays intermedios son de longitud 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);
}
Encontrar el punto del segmento más próximo:
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;
}
Este es un método simple de prueba si el punto T, que sabemos que está en la misma gran círculo como A y B, está en el segmento más corto de este gran círculo. Sin embargo, hay métodos más eficientes de hacerlo:
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;
}
Otros consejos
Trate distancia de un punto a un círculo mayor , desde ask Dr. . Matemáticas. Todavía es necesario transformar la longitud / latitud a coordenadas esféricas y la escala para el radio de la tierra, pero esto parece una buena dirección.
Este es el código completo para la respuesta aceptada como Ideone violín (que se encuentra aquí ):
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();
}
}
Funciona bien para los datos comentado:
//line
double [] a = {50.174315,19.054743};
double [] b = {50.176019,19.065042};
//point
double [] c = {50.184373,19.054657};
nodo más cercano es: 50.17493121381319,19.05846668493702
Pero tengo un problema con estos datos:
double [] a = {52.00118, 17.53933};
double [] b = {52.00278, 17.54008};
//point
double [] c = {52.008308, 17.542927};
nodo más cercano es: 52.00834987257176,17.542691313436357 que está mal
.Creo que de línea especificado por dos puntos no es un segmento cerrado.
Si alguien necesita que esta es la respuesta loleksy portado a 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;
}
Para una distancia de hasta unos pocos miles de metros que simplificaría el problema de esfera en avión. Entonces, la cuestión es bastante simple como un cálculo fácil triángulo se puede utilizar:
Tenemos los puntos A y B y buscamos una distancia X a la línea AB. A continuación:
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;
La distancia más corta entre dos puntos de una esfera es la parte más pequeña de la gran círculo que pasa por los dos puntos. Estoy seguro de que ya sabe esto. Hay una pregunta similar aquí http://www.physicsforums.com/ archive / index.php / t-178252.html que puede ayudar a modelar se mathmatically.
No estoy seguro de cómo será la probabilidad de obtener un código de ejemplo de este, para ser honesto.
Estoy buscando básicamente lo mismo en este momento, excepto que en sentido estricto no se preocupan por tener un segmento de un círculo, sino que sólo quieren la distancia a cualquier punto de la circunferencia completa.
Dos enlaces Actualmente estoy investigando:
Esta página menciona "Cross-track distancia", que básicamente parece ser lo que busca.
Además, en el siguiente hilo en la lista de correo PostGIS, el intento parece (1) determinar el punto más cercano en el gran círculo con la misma fórmula usada para la línea de distancia en un 2D-plano (con line_locate_point PostGIS') y, a continuación (2) cálculo de la distancia entre esa y el tercer punto en un esferoide. No tengo idea de si el paso matemáticamente (1) es correcta, pero me sorprendería.
http://postgis.refractions.net/pipermail /postgis-users/2009-July/023903.html
Por último, acabo de ver que los siguientes ligados en "relacionada":
distancia desde el punto a la línea de gran función círculo no funciona bien.
