Como calcular a distância de um segmento de ponto a um linha, em uma esfera?
https://stackoverflow.com/questions/1299567
-
18-09-2019 - |
italiano
english
français
española
中国
日本の
العربية
Deutsch
한국어
Português
RussianPergunta
Eu tenho um segmento de linha (grande parte do círculo) na Terra. O segmento de linha é definido pelas coordenadas de suas extremidades. Obviamente, dois pontos definem dois segmentos de linha, então assuma que estou interessado no mais curto.
Recebi um terceiro ponto e estou procurando a distância (menor) entre a linha e o ponto.
Todas as coordenadas são dadas em longitude latitude (WGS 84).
Como faço para calcular a distância?
Uma solução em qualquer linguagem de programação razoável serve.
Solução
Aqui está minha própria solução, com base na ideia em Pergunte ao Dr. Math. Eu ficaria feliz em ver seu feedback.
Isenção de responsabilidade primeiro. Esta solução está correta para as esferas. A Terra não é uma esfera, e o sistema de coordenadas (WGS 84) não assume que é uma esfera. Portanto, isso é apenas uma aproximação, e eu realmente não posso estimar é um erro. Além disso, para distâncias muito pequenas, provavelmente também é possível obter uma boa aproximação assumindo que tudo é apenas um coplanar. Mais uma vez, não sei como as distâncias "pequenas" devem ser.
Agora para os negócios. Vou chamar as extremidades das linhas A, B e o terceiro ponto C. Basicamente, o algoritmo é: para:
- Converta as coordenadas primeiro em coordenadas cartesianas (com a origem no centro da terra) - por exemplo, aqui.
Calcule T, o ponto da linha AB que está mais próximo de C, usando os três produtos de vetores a seguir:
G = a x b
F = c x g
T = g x f
Normalize T e multiplique pelo raio da terra.
- Converta T de volta à longitude Latitude.
- Calcule a distância entre t e c - por exemplo, aqui.
Essas etapas são suficientes se você estiver procurando a distância entre C e o Grande Círculo Definado por A e B. Se como eu, você está interessado na distância entre C e o segmento de linha mais curto, você precisa dar o passo extra para verificar que T está realmente nesse segmento. Se não for, então necessariamente o ponto mais próximo é uma das extremidades A ou B - a maneira mais fácil é verificar qual.
Em termos gerais, a idéia por trás dos três produtos vetoriais é a seguinte. O primeiro (g) nos dá o plano do grande círculo de A e B (para que o plano contendo A, B e a Origem). O segundo (f) nos dá o grande círculo que passa por C e é perpendicular ao G. Então T é a interseção dos grandes círculos definidos por F e G, trazidos ao comprimento correto por normalização e multiplicação por R.
Aqui está algum código Java parcial para fazê -lo.
Encontrando o ponto mais próximo no grande círculo. As entradas e a saída são matrizes de comprimento-2. As matrizes intermediárias são de comprimento 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);
}
Encontrando o ponto mais próximo no 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 é um método simples de teste se o ponto T, que sabemos que está no mesmo grande círculo que A e B, está no segmento mais curto deste grande círculo. No entanto, existem métodos mais eficientes para fazê -lo:
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;
}
Outras dicas
Tentar Distância de um ponto a um grande círculo, de Ask Dr. Math. Você ainda precisa transformar a longitude/latitude em coordenadas esféricas e escalar para o raio da Terra, mas isso parece uma boa direção.
Este é o código completo para a resposta aceita como violino ideal (encontrado aqui):
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 bem para dados comentados:
//line
double [] a = {50.174315,19.054743};
double [] b = {50.176019,19.065042};
//point
double [] c = {50.184373,19.054657};
O nó mais próximo é: 50.17493121381319.19.05846668493702
Mas tenho problemas com esses dados:
double [] a = {52.00118, 17.53933};
double [] b = {52.00278, 17.54008};
//point
double [] c = {52.008308, 17.542927};
O nó mais próximo é: 52.00834987257176.17.542691313436357, o que está errado.
Eu acho que essa linha especificada por dois pontos não é um segmento fechado.
Se alguém precisar, essa é a resposta loleksy portada para 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 uma distância de até alguns milhares de metros, simplificaria o problema de esfera em avião. Então, o problema é simplesmente como um cálculo fácil do triângulo pode ser usado:
Temos pontos A e B e procuram uma distância x para alinhar AB. Então:
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;
A distância mais curta entre dois pontos em uma esfera é o lado menor do grande círculo que passa pelos dois pontos. Tenho certeza que você já sabe disso. Há uma pergunta semelhante aqui http://www.physicsforums.com/archive/index.php/t-178252.html Isso pode ajudá -lo a modelá -lo matematicamente.
Não tenho certeza de quão provável você tem um exemplo codificado disso, para ser honesto.
Estou basicamente procurando a mesma coisa agora, exceto que estritamente não me importo em ter um segmento de um grande círculo, mas só quero a distância de qualquer ponto no círculo completo.
Dois links que estou investigando atualmente:
Esta página menciona a "distância cruzada", o que basicamente parece ser o que você está procurando.
Além disso, no encadeamento a seguir na lista de discussão do PostGIS, a tentativa parece (1) determinar o ponto mais próximo do grande círculo com a mesma fórmula usada para distância de linha em um plano 2D (com pós-line_locate_point) e depois (2) calcular a distância entre isso e o terceiro ponto em um esferóide. Não faço ideia se a etapa matematicamente (1) está correta, mas ficaria surpreso.
http://postgis.refrações.net/pipermail/postgis-users/2009-july/023903.html
Finalmente, acabei de ver que o seguinte vinculado em "Related":
Distância de ponto a linha Great Circle Função não está funcionando corretamente.
