language agnostic - terrestres - ¿Cómo calcular la distancia de un punto a un segmento de línea, en una esfera?
formula para calcular distancia entre dos puntos geograficos (7)
Tengo un segmento de línea (gran círculo) en la tierra. El segmento de línea está definido por las coordenadas de sus extremos. Obviamente, dos puntos definen dos segmentos de línea, así que supongamos que estoy interesado en el más corto.
Me dan un tercer punto, y estoy buscando la distancia (más corta) entre la línea y el punto.
Todas las coordenadas se dan en longitud / latitud (WGS 84).
¿Cómo calculo la distancia?
Una solución en cualquier lenguaje de programación razonable servirá.
Aquí está mi propia solución, basada en la idea de preguntar al Dr. Math . Estaría feliz de ver tus comentarios.
Descargo de responsabilidad primero. Esta solución es correcta para las esferas. La Tierra no es una esfera, y el sistema de coordenadas (WGS 84) no supone que sea una esfera. Entonces esto es solo una aproximación, y realmente no puedo estimar el error. Además, para distancias muy pequeñas, probablemente también sea posible obtener una buena aproximación asumiendo que todo es simplemente coplanaire. De nuevo, no sé cuán "pequeñas" deben ser las distancias.
Ahora a los negocios. Voy a llamar a los extremos de las líneas A, B y el tercer punto C. Básicamente, el algoritmo es:
- convierta las coordenadas primero en coordenadas cartesianas (con el origen en el centro de la tierra), por ejemplo, aquí .
Calcule T, el punto en la línea AB que está más cerca de C, usando los siguientes 3 productos vectoriales:
G = A x B
F = C x G
T = G x F
Normaliza T y multiplica por el radio de la tierra.
- Convierta T de nuevo a longitud / latitud.
- Calcule la distancia entre T y C, por ejemplo aquí .
Estos pasos son suficientes si está buscando la distancia entre C y el gran círculo definido por A y B. Si, como a mí, le interesa la distancia entre C y el segmento de línea más corto, debe tomar el paso adicional de verificar que T está de hecho en este segmento. Si no lo es, entonces necesariamente el punto más cercano es uno de los extremos A o B - la manera más fácil es verificar cuál.
En términos generales, la idea detrás de los tres productos vectoriales es la siguiente. El primero (G) nos da el plano del gran círculo de A y B (por lo tanto, el plano que contiene A, B y el origen). El segundo (F) nos da el gran círculo que pasa por C y es perpendicular al G. Entonces T es la intersección de los grandes círculos definidos por F y G, llevados a la longitud correcta por normalización y multiplicación por R.
Aquí hay un código parcial de Java para hacerlo.
Encontrar el punto más cercano en el gran círculo. Las entradas y salidas son matrices de longitud 2. Los arreglos 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 más cercano en el segmento:
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 el mismo círculo grande que A y B, está en el segmento más corto de este gran círculo. Sin embargo, existen métodos más eficientes para 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;
}
Básicamente estoy buscando lo mismo en este momento, excepto que estrictamente hablando no me importa tener un segmento de un gran círculo, sino que simplemente quiero la distancia a cualquier punto del círculo completo.
Dos enlaces que estoy investigando actualmente:
Esta página menciona la "distancia entre pistas", que básicamente parece ser lo que estás buscando.
Además, en el siguiente hilo de 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 utilizada para la distancia de línea en un plano 2D (con el punto de ubicación de línea de PostGIS), y luego (2) calcular la distancia entre ese y el tercer punto en un esferoide. No tengo idea si matemáticamente el paso (1) es correcto, pero me sorprendería.
http://postgis.refractions.net/pipermail/postgis-users/2009-July/023903.html
Finalmente, acabo de ver que lo siguiente se vincula bajo "Relacionados":
La distancia desde la función de gran círculo de punto a línea no funciona bien.
Este es el código completo para la respuesta aceptada como ideone fiddle (se encuentra 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://.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://.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://.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 comentados:
//line
double [] a = {50.174315,19.054743};
double [] b = {50.176019,19.065042};
//point
double [] c = {50.184373,19.054657};
El nodo más cercano es: 50.17493121381319,19.05846668493702
Pero tengo un problema con esta información:
double [] a = {52.00118, 17.53933};
double [] b = {52.00278, 17.54008};
//point
double [] c = {52.008308, 17.542927};
El nodo más cercano es: 52.00834987257176,17.542691313436357 que es incorrecto.
Creo que la línea especificada por dos puntos no es un segmento cerrado.
La distancia más corta entre dos puntos en una esfera es el lado más pequeño del gran círculo que pasa por los dos puntos. Estoy seguro de que ya lo sabes. Hay una pregunta similar aquí http://www.physicsforums.com/archive/index.php/t-178252.html que puede ayudarlo a modelar matemáticamente.
No estoy seguro de qué tan probable sea obtener un ejemplo codificado de esto, para ser honesto.
Para una distancia de hasta unos pocos miles de metros, simplificaría el problema de esfera a plano. Entonces, el problema es bastante simple ya que se puede usar un cálculo de triángulo fácil:
Tenemos los puntos A y B y buscamos una distancia X a la línea AB. Entonces:
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;
Pruebe la distancia desde un punto a un gran círculo , desde Ask Dr. Math. Aún necesita transformar la longitud / latitud a coordenadas esféricas y escalar para el radio de la Tierra, pero parece una buena dirección.
Si alguien lo necesita, esta es la respuesta loleksportada 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://.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://.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://.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;
}