Trilateração usando 3 pontos de latitude / longitude e 3 distâncias?

34

Quero descobrir um local de destino desconhecido (coordenadas de latitude e longitude). Existem 3 pontos conhecidos (pares de coordenadas de latitude e longitude) e, para cada ponto, uma distância em quilômetros até o local de destino. Como posso calcular as coordenadas do local de destino?

Por exemplo, digamos que tenho os seguintes pontos de dados

37.418436,-121.963477   0.265710701754km
37.417243,-121.961889   0.234592423446km
37.418692,-121.960194   0.0548954278262km

O que eu gostaria é o que é a matemática de uma função que aceita isso como entrada e retorna 37.417959, -121.961954 como saída.

Entendo como calcular a distância entre dois pontos, de http://www.movable-type.co.uk/scripts/latlong.html Entendo o princípio geral de que, com três círculos como esses, você obtém exatamente um ponto de sobreposição. O que eu estou confuso é a matemática necessária para calcular esse ponto com essa entrada.

nohat
fonte
Aqui está uma página que mostra a matemática para encontrar o centro das três coordenadas. Talvez possa ajudar de alguma forma. < mathforum.org/library/drmath/view/68373.html >
Jon Bringhurst
11
Isso precisa estar na esfera / esferóide ou um algoritmo planar está bom?
fmark 23/07/10
11
Não posso lhe dar a resposta, mas acho que posso apontá-lo na direção certa. Três coordenadas = três pontos centrais. Três distâncias = três círculos. Dois círculos que se cruzam, podem ter a possibilidade de nenhuma / uma / duas soluções. Três círculos podem ter nenhuma / uma / ou uma área como solução. Obtenha a fórmula do círculo para os três círculos e resolva-a com Sistemas de Equações / Álgebra.
23410 CrazyEnigma
Na verdade, você nem precisa de sistemas para resolver este. Existem uma ou duas possibilidades, mas como você tem um valor de distância, pode separar a resposta correta.
George Silva
11
+1 Esta é uma boa pergunta. No começo, achei que uma solução poderia ser facilmente encontrada com o google, mas aparentemente não. Talvez o problema possa ser mais generalizado: dados N pontos com cada ponto tendo não apenas uma distância, mas também uma margem de erro, encontre a elipse de confiança.
22410 Kirk Kuykendall

Respostas:

34

Depois de analisar a Wikipedia e a mesma pergunta / resposta no StackOverflow , imaginei que iria dar uma facada nele e tentar preencher as lacunas.

Primeiro, não sei onde você obteve a saída, mas ela parece estar errada. Plotei os pontos no ArcMap, armazenei-os em buffer nas distâncias especificadas, executei a interseção nos buffers e depois capturei o vértice da interseção para obter as soluções. Sua saída proposta é o ponto em verde. Calculei o valor na caixa de texto explicativo, que é cerca de 3 metros do que o ArcMap deu para a solução derivada da interseção.

texto alternativo

A matemática na página da Wikipédia não é muito ruim, basta ocultar suas coordenadas geodésicas ao ECEF cartesiano, que pode ser encontrado aqui . os termos a / x + h podem ser substituídos pelo raio da esfera autálica, se você não estiver usando um elipsóide.

Provavelmente mais fácil, basta fornecer um código bem documentado (?), Então aqui está em python

import math
import numpy

#assuming elevation = 0
earthR = 6371
LatA = 37.418436
LonA = -121.963477
DistA = 0.265710701754
LatB = 37.417243
LonB = -121.961889
DistB = 0.234592423446
LatC = 37.418692
LonC = -121.960194
DistC = 0.0548954278262

#using authalic sphere
#if using an ellipsoid this step is slightly different
#Convert geodetic Lat/Long to ECEF xyz
#   1. Convert Lat/Long to radians
#   2. Convert Lat/Long(radians) to ECEF
xA = earthR *(math.cos(math.radians(LatA)) * math.cos(math.radians(LonA)))
yA = earthR *(math.cos(math.radians(LatA)) * math.sin(math.radians(LonA)))
zA = earthR *(math.sin(math.radians(LatA)))

xB = earthR *(math.cos(math.radians(LatB)) * math.cos(math.radians(LonB)))
yB = earthR *(math.cos(math.radians(LatB)) * math.sin(math.radians(LonB)))
zB = earthR *(math.sin(math.radians(LatB)))

xC = earthR *(math.cos(math.radians(LatC)) * math.cos(math.radians(LonC)))
yC = earthR *(math.cos(math.radians(LatC)) * math.sin(math.radians(LonC)))
zC = earthR *(math.sin(math.radians(LatC)))

P1 = numpy.array([xA, yA, zA])
P2 = numpy.array([xB, yB, zB])
P3 = numpy.array([xC, yC, zC])

#from wikipedia
#transform to get circle 1 at origin
#transform to get circle 2 on x axis
ex = (P2 - P1)/(numpy.linalg.norm(P2 - P1))
i = numpy.dot(ex, P3 - P1)
ey = (P3 - P1 - i*ex)/(numpy.linalg.norm(P3 - P1 - i*ex))
ez = numpy.cross(ex,ey)
d = numpy.linalg.norm(P2 - P1)
j = numpy.dot(ey, P3 - P1)

#from wikipedia
#plug and chug using above values
x = (pow(DistA,2) - pow(DistB,2) + pow(d,2))/(2*d)
y = ((pow(DistA,2) - pow(DistC,2) + pow(i,2) + pow(j,2))/(2*j)) - ((i/j)*x)

# only one case shown here
z = numpy.sqrt(pow(DistA,2) - pow(x,2) - pow(y,2))

#triPt is an array with ECEF x,y,z of trilateration point
triPt = P1 + x*ex + y*ey + z*ez

#convert back to lat/long from ECEF
#convert to degrees
lat = math.degrees(math.asin(triPt[2] / earthR))
lon = math.degrees(math.atan2(triPt[1],triPt[0]))

print lat, lon
wwnick
fonte
11
Eu estava indo para montar uma resposta semelhante, mas agora não há necessidade! Recebe meu voto positivo.
Wrass
entorpecido para o resgate! Compila quando 'triPt' é substituído por 'triLatPt', mas, caso contrário, retorna 37.4191023738 -121.960579208. Bom trabalho #
WolfOdrade 28/07
Bom trabalho! Se eu substituir o sistema de coordenadas geográficas por um sistema de coordenadas [cartesiano] local, isso ainda funcionará?
Zengr 03/10
para aqueles no domínio c ++ .. combinaram um bem rápido pastebin.com/9Dur6RAP
raaj
2
Obrigado @wwnick! Portamos isso para JavaScript (destinado ao Node, mas pode ser facilmente convertido para funcionar no navegador). gist.github.com/dav-/bb7103008cdf9359887f
DC_
6

Não tenho certeza se estou sendo ingênuo, mas se você colocar cada ponto em buffer por tamanho e, em seguida, cruzar todos os três círculos que o levariam ao local correto?

Você pode calcular a interseção usando APIs espaciais. Exemplos:

  • GeoScript
  • Conjunto de Topologia Java
  • NET Topology Suite
  • GEOS
George Silva
fonte
11
Exatamente, ele está interessado nas fórmulas para obter esse ponto de interseção.
Vinko Vrsalovic 22/07/10
Usando uma API espacial, você pode fazer isso sem usar matemática pura.
George Silva
11
@ George Você pode dar um exemplo dessa API?
nohat 23/07
Postagem editada para refletir a solicitação de nohat.
George Silva
+1, um bom pensamento lateral, mesmo que não seja o mais eficiente em termos de computação!
fmark
2

As notas a seguir usam geometria planarítmica (ou seja, você teria que projetar suas coordenadas em um sistema de coordenadas local apropriado).

Meu raciocínio, com um exemplo trabalhado em Python, segue:

Pegue 2 dos pontos de dados (chame-os ae b). Ligue para o nosso ponto alvo x. Já sabemos as distâncias axe bx. Podemos calcular a distância abusando o teorema de Pitágoras.

>>> import math
>>> a = (1, 4)
>>> b = (3, 6)
>>> dist_ax = 3
>>> dist_bx = 5.385
# Pythagoras's theorem
>>> dist_ab = math.sqrt(abs(a[0]-b[0])**2 + abs(a[1]-b[1])**2)
>>> dist_ab
2.8284271247461903

Agora, você pode trabalhar os ângulos dessas linhas:

>>> angle_abx = math.acos((dist_bx * dist_bx + dist_ab * dist_ab - dist_ax * dist_ax)/(2 * dist_bx * dist_ab))
>>> math.degrees(angle_abx)
23.202973815040256
>>> angle_bax = math.acos((dist_ax * dist_ax + dist_ab * dist_ab - dist_bx * dist_bx)/(2 * dist_ax * dist_ab))
>>> math.degrees(angle_bax)
134.9915256259537
>>> angle_axb = math.acos((dist_ax * dist_ax + dist_bx * dist_bx - dist_ab * dist_ab)/(2 * dist_ax * dist_bx))
>>> math.degrees(angle_axb)
21.805500559006095

Infelizmente, estou com pouco tempo para responder a você. No entanto, agora que você conhece os ângulos, pode calcular dois locais possíveis x. Em seguida, usando o terceiro ponto c, você pode calcular qual local está correto.

fmark
fonte
2

Isso pode funcionar. Rapidamente novamente em python, você pode colocar isso no corpo de uma função xN, yN = coordenadas dos pontos, r1 & r2 = valores do raio

dX = x2 - x1
dY = y2 - y1

centroidDistance = math.sqrt(math.pow(e,2) + math.pow(dY,2)) #distance from centroids
distancePL = (math.pow(centroidDistance,2) + (math.pow(r1,2) - math.pow(r2,2))) / (2 * centroidDistance) #distance from point to a line splitting the two centroids

rx1 = x1 + (dX *k)/centroidDistance + (dY/centroidDistance) * math.sqrt(math.pow(r1,2) - math.pow(distancePL,2))
ry1 = y1 + (dY*k)/centroidDistance - (dX /centroidDistance) * math.sqrt(math.pow(r1,2) - math.pow(distancePL,2))

rx2 = x1 + (dX *k)/centroidDistance - (dY/centroidDistance) * math.sqrt(math.pow(r1,2) - math.pow(distancePL,2))
ry2 = y1 + (dY*k)/centroidDistance + (dX /centroidDistance) * math.sqrt(math.pow(r1,2) - math.pow(distancePL,2))

Os valores rx e ry são os valores de retorno (devem estar em uma matriz) dos dois pontos de interseção em um círculo, se isso ajudar a esclarecer as coisas.

Faça isso nos 2 primeiros círculos e depois no primeiro e no último. Se um dos resultados da primeira iteração se comparar com os resultados da segunda (provavelmente com alguma tolerância), você terá o ponto de interseção. Não é uma ótima solução, especialmente quando você começa a adicionar mais do que pontos no processo, mas é a mais simples que posso ver sem resolver um sistema de equações.

WolfOdrade
fonte
O que são 'e' e 'k' no seu código?
ReinierDG 14/09
Não me lembro :-) A resposta de wwnick é mais parecida com algo que você gostaria de implementar se tivesse apenas três círculos.
WolfOdrade
1

Você pode usar a API espacial do postgis (funções St_Intersection, St_buffer). Como fmark observou, você também deve se lembrar que o Postgis usa algoritmos planares, mas para pequenas áreas, o uso de projeção equidistante não gera muitos erros.

stachu
fonte
O PostGIS pode fazer cálculos esferoidais usando o GEOGRAPHYtipo e não o GEOMETRYtipo.
fmark 25/07/10
1

Faça isso na linguagem PHP:

// assumindo elevação = 0
$ earthR = 6371; // em km (= 3959 em milhas)

$ LatA = 37,418436;
$ LonA = -121,963477;
$ DistA = 0,265710701754;

$ LatB = 37,417243;
$ LonB = -121,961889;
$ DistB = 0,234592423446;

$ LatC = 37,418692;
$ LonC = -121.960194;
$ DistC = 0,0548954278262;

/ *
#using esfera authalic
#se usando um elipsóide, esta etapa é um pouco diferente
#Converte geodésico Lat / Long em ECEF xyz
# 1. Converter lat / long para radianos
# 2. Converter lat / long (radianos) para ECEF
* /
$ xA = $ earthR * (cos (deg2rad ($ LatA)) * cos (deg2rad ($ LonA)));
$ yA = $ earthR * (cos (deg2rad ($ LatA)) * sin (deg2rad ($ LonA)));
$ zA = $ earthR * (sin (deg2rad ($ LatA)));

$ xB = $ earthR * (cos (deg2rad ($ LatB)) * cos (deg2rad ($ LonB)));
$ yB = $ earthR * (cos (deg2rad ($ LatB)) * sin (deg2rad ($ LonB)));
$ zB = $ earthR * (sin (deg2rad ($ LatB)));

$ xC = $ earthR * (cos (deg2rad ($ LatC)) * cos (deg2rad ($ LonC)));
$ yC = $ earthR * (cos (deg2rad ($ LatC)) * sin (deg2rad ($ LonC)));
$ zC = $ earthR * (sin (deg2rad ($ LatC)));

/ *
INSTALAR:
Instalação do sudo pear Math_Vector-0.7.0
Instalação do sudo pear Math_Matrix-0.8.7
* /
// Inclui PEAR :: Math_Matrix
// /usr/share/php/Math/Matrix.php
// include_path = ".: / usr / local / php / pear /"
require_once 'Math / Matrix.php';
require_once 'Math / Vector.php';
require_once 'Math / Vector3.php';


$ P1vector = new Math_Vector3 (matriz ($ xA, $ yA, $ zA));
$ P2vector = novo Math_Vector3 (matriz ($ xB, $ yB, $ zB));
$ P3vector = new Math_Vector3 (matriz ($ xC, $ yC, $ zC));

#from wikipedia: http://en.wikipedia.org/wiki/Trilateration
#transform para obter o círculo 1 na origem
#transform para obter o círculo 2 no eixo x

// CALC EX
$ P2minusP1 = Math_VectorOp :: substract ($ P2vector, $ P1vector);
$ l = novo Math_Vector ($ P2minusP1);
$ P2minusP1_length = $ l-> length ();
$ norm = new Math_Vector3 (matriz ($ P2minusP1_length, $ P2minusP1_length, $ P2minusP1_length));
$ d = $ norma; // salva calc D
$ ex = Math_VectorOp :: divide ($ P2minusP1, $ norm);
// echo "ex:". $ ex-> toString (). "\ n";
$ ex_x = floatval ($ ex -> _ tupla-> getData () [0]);
$ ex_y = floatval ($ ex -> _ tupla-> getData () [1]);
$ ex_z = floatval ($ ex -> _ tupla-> getData () [2]);
$ ex = novo Math_Vector3 (matriz ($ ex_x, $ ex_y, $ ex_z));

// CALC i
$ P3minusP1 = Math_VectorOp :: substract ($ P3vector, $ P1vector);
$ P3minusP1_x = floatval ($ P3minusP1 -> _ tupla-> getData () [0]);
$ P3minusP1_y = floatval ($ P3minusP1 -> _ tupla-> getData () [1]);
$ P3minusP1_z = floatval ($ P3minusP1 -> _ tupla-> getData () [2]);
$ P3minusP1 = novo Math_Vector3 (matriz ($ P3minusP1_x, $ P3minusP1_y, $ P3minusP1_z));
$ i = Math_VectorOp :: dotProduct ($ ex, $ P3minusP1);
// echo "i = $ i \ n";

// CALC EY
$ iex = Math_VectorOp :: scale ($ i, $ ex);
// echo "iex =". $ iex-> toString (). "\ n";
$ P3P1iex = Math_VectorOp :: substract ($ P3minusP1, $ iex);
// echo "P3P1iex =". $ P3P1iex-> toString (). "\ n";
$ l = novo Math_Vector ($ P3P1iex);
$ P3P1iex_length = $ l-> length ();
$ norm = new Math_Vector3 (array ($ P3P1iex_length, $ P3P1iex_length, $ P3P1iex_length));
// echo "norm:". $ norm-> toString (). "\ n";
$ ey = Math_VectorOp :: divide ($ P3P1iex, $ norm);
// echo "ey =". $ ey-> toString (). "\ n";
$ ey_x = floatval ($ ey -> _ tupla-> getData () [0]);
$ ey_y = floatval ($ ey -> _ tupla-> getData () [1]);
$ ey_z = floatval ($ ey -> _ tupla-> getData () [2]);
$ ey = new Math_Vector3 (matriz ($ ey_x, $ ey_y, $ ey_z));

// CALC EZ
$ ez = Math_VectorOp :: crossProduct ($ ex, $ ey);
// echo "ez =". $ ez-> toString (). "\ n";

// CALC D
// faça isso antes
$ d = valor flutuante ($ d -> _ tupla-> getData () [0]);
// echo "d = $ d \ n";

// CALC J
$ j = Math_VectorOp :: dotProduct ($ ey, $ P3minusP1);
// echo "j = $ j \ n";

#from wikipedia
#plug and chug usando os valores acima
$ x = (pow ($ DistA, 2) - pow ($ DistB, 2) + pow ($ d, 2)) / (2 * $ d);
$ y = ((pow ($ DistA, 2) - pow ($ DistC, 2) + pow ($ i, 2) + pow ($ j, 2)) / (2 * $ j)) - (($ i / $ j) * $ x);

# Apenas um caso mostrado aqui
$ z = sqrt (pow ($ DistA, 2) - pow ($ x, 2) - pow ($ y, 2));

// echo "x = $ x - y = $ y - z = $ z \ n";

#triPt é uma matriz com ECEF x, y, z do ponto de trilateração
$ xex = Math_VectorOp :: scale ($ x, $ ex);
$ yey = Math_VectorOp :: scale ($ y, $ ey);
$ zez = Math_VectorOp :: scale ($ z, $ ez);

// CALC $ triPt = $ P1vector + $ xex + $ yey + $ zez;
$ triPt = Math_VectorOp :: add ($ P1vector, $ xex);
$ triPt = Math_VectorOp :: add ($ triPt, $ yey);
$ triPt = Math_VectorOp :: add ($ triPt, $ zez);
// echo "triPt =". $ triPt-> toString (). "\ n";
$ triPt_x = floatval ($ triPt -> _ tupla-> getData () [0]);
$ triPt_y = floatval ($ triPt -> _ tupla-> getData () [1]);
$ triPt_z = floatval ($ triPt -> _ tupla-> getData () [2]);


#converter de volta para lat / long do ECEF
#convert em graus
$ lat = rad2deg (asin ($ triPt_z / $ earthR));
$ lon = rad2deg (atan2 ($ triPt_y, $ triPt_x));

echo $ lat. ','. $ lon;
fquinto
fonte