Trilateracja przy użyciu 3 punktów szerokości / długości geograficznej i 3 odległości?

Chcę znaleźć nieznaną lokalizację docelową (współrzędne szerokości i długości geograficznej). Istnieją 3 znane punkty (pary współrzędnych szerokości i długości geograficznej) i dla każdego punktu odległość w kilometrach od miejsca docelowego. Jak mogę obliczyć współrzędne docelowej lokalizacji?

Powiedzmy na przykład, że mam następujące punkty danych

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

To, co chciałbym, to matematyka dla funkcji, która przyjmuje to jako dane wejściowe i zwraca 37,417959, -121,961954 jako dane wyjściowe.

Rozumiem, jak obliczyć odległość między dwoma punktami, na stronie http://www.movable-type.co.uk/scripts/latlong.html Rozumiem ogólną zasadę, że przy takich trzech okręgach otrzymujesz dokładnie jeden punkt nakładania się. Mętnie rozumiem matematykę potrzebną do obliczenia tego punktu za pomocą tych danych wejściowych.

Po tym, jak rozejrzałem się po Wikipedii i tym samym pytaniu / odpowiedzi na StackOverflow , pomyślałem, że zrobię sobie dźgnięcie i spróbuję uzupełnić luki.

Po pierwsze, Nie jestem pewien, skąd masz wynik, ale wydaje się, że jest nieprawidłowy. Wykreśliłem punkty w ArcMap, zbuforowałem je do określonych odległości, pobiegłem przeciąć na buforach, a następnie uchwyciłem wierzchołek przecięcia, aby uzyskać rozwiązania. Proponowany wynik jest zielony. Obliczyłem wartość w oknie objaśnień, czyli około 3 metry tego, co ArcMap dał dla rozwiązania uzyskanego z przecięcia.

Matematyka na stronie wikipedii nie jest taka zła, wystarczy zakryć swoje współrzędne geodezyjne kartezjańskim ECEF, który można znaleźć tutaj . warunki a / x + h można zastąpić promieniem kuli authalicznej, jeśli nie używasz elipsoidy.

Prawdopodobnie najłatwiej jest po prostu podać ci dobrze (?) Udokumentowany kod, więc tutaj jest w Pythonie

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

Nie jestem pewien, czy jestem naiwny, ale czy buforujesz każdy punkt według rozmiaru, a następnie przecinasz wszystkie trzy koła, które zapewniłyby ci poprawną lokalizację?

Możesz obliczyć skrzyżowanie za pomocą przestrzennych interfejsów API. Przykłady:

• GeoScript
• Java Topology Suite
• NET Topology Suite
• GEOS

Poniższe uwagi wykorzystują geometrię planaryczną (tj. Musiałbyś rzutować swoje współrzędne na odpowiedni lokalny układ współrzędnych).

Moje rozumowanie z działającym przykładem w Pythonie wygląda następująco:

Weź 2 punkty danych (zadzwoń do nich ai b). Zadzwoń do naszego punktu docelowego x. Znamy już odległości axi bx. Odległość możemy obliczyć za abpomocą twierdzenia Pitagorasa.

>>> 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

Teraz możesz obliczyć kąty tych linii:

>>> 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

Niestety brakuje mi czasu na udzielenie odpowiedzi, jednak teraz znasz kąty, możesz obliczyć dwie możliwe lokalizacje x. Następnie, korzystając z trzeciego punktu c, możesz obliczyć, która lokalizacja jest poprawna.

To może zadziałać. Szybko ponownie w pythonie, możesz umieścić to w ciele funkcji xN, yN = współrzędne punktów, r1 i r2 = wartości promienia

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))

Wartości rx i ry to wartości zwracane (powinny znajdować się w tablicy) dwóch punktów przecięcia na okręgu, jeśli to pomaga wyjaśnić.

Zrób to dla pierwszych 2 kręgów, a następnie dla pierwszego i ostatniego. Jeśli którykolwiek z wyników z pierwszej iteracji porówna się z wynikami z drugiego (prawdopodobnie w ramach pewnej tolerancji), to masz punkt przecięcia. Nie jest to świetne rozwiązanie, zwłaszcza gdy zaczynasz dodawać do procesu więcej niż punkty, ale jest to najprostsze, jakie widzę bez rozwiązywania układu równań.

Możesz używać przestrzennego API z postgis (funkcje St_Intersection, St_buffer). Jak zauważył fmark, należy również pamiętać, że Postgis stosuje algorytmy planarne, ale w przypadku małych obszarów stosowanie równomiernego przesuwania nie wprowadza dużego błędu.

Zrób to w języku PHP:

// przy założeniu rzędnej = 0
$ earthR = 6371; // w km (= 3959 w milach)

$ 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;

/ *
#z wykorzystaniem kuli autologicznej
#if przy użyciu elipsoidy ten krok jest nieco inny
# Przelicz geodezyjny Lat / Long na ECEF xyz
# 1. Konwertuj Lat / Long na radiany
# 2. Konwertuj Lat / Long (radiany) na 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)));

/ *
ZAINSTALOWAĆ:
sudo gruszka zainstalować Math_Vector-0.7.0
sudo gruszka zainstalować Math_Matrix-0.8.7
* /
// Dołącz PEAR :: Math_Matrix
// /usr/share/php/Math/Matrix.php
// include_path = ".: / usr / local / php / pear /"
requ_once 'Math / Matrix.php';
requ_once 'Math / Vector.php';
requ_once 'Math / Vector3.php';


$ P1vector = nowy Math_Vector3 (tablica ($ xA, $ yA, $ zA));
$ P2vector = nowy Math_Vector3 (tablica ($ xB, $ yB, $ zB));
$ P3vector = nowy Math_Vector3 (tablica ($ xC, $ yC, $ zC));

# z wikipedii: http://en.wikipedia.org/wiki/Trilateration
#transformuj, aby uzyskać krąg 1 w punkcie początkowym
#transformuj, aby uzyskać okrąg 2 na osi x

// CALC EX
$ P2minusP1 = Math_VectorOp :: subtract ($ P2vector, $ P1vector);
$ l = nowy Math_Vector ($ P2minusP1);
$ P2minusP1_length = $ l-> length ();
$ norm = nowy Math_Vector3 (tablica ($ P2minusP1_length, $ P2minusP1_length, $ P2minusP1_length));
$ d = $ norma; // zapisz oblicz D
$ ex = Math_VectorOp :: divide ($ P2minusP1, $ norm);
// echo "ex:". $ ex-> toString (). "\ n";
$ ex_x = floatval ($ ex -> _ tuple-> getData () [0]);
$ ex_y = floatval ($ ex -> _ tuple-> getData () [1]);
$ ex_z = floatval ($ ex -> _ tuple-> getData () [2]);
$ ex = nowy Math_Vector3 (tablica ($ ex_x, $ ex_y, $ ex_z));

// CALC i
$ P3minusP1 = Math_VectorOp :: subtract ($ P3vector, $ P1vector);
$ P3minusP1_x = zmiennoprzecinkowy ($ P3minusP1 -> _ tuple-> getData () [0]);
$ P3minusP1_y = floatval ($ P3minusP1 -> _ tuple-> getData () [1]);
$ P3minusP1_z = floatval ($ P3minusP1 -> _ tuple-> getData () [2]);
$ P3minusP1 = nowy Math_Vector3 (tablica ($ 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 :: subtract ($ P3minusP1, $ iex);
// echo "P3P1iex =". $ P3P1iex-> toString (). "\ n";
$ l = nowy Math_Vector ($ P3P1iex);
$ P3P1iex_length = $ l-> length ();
$ norm = nowy Math_Vector3 (tablica ($ 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 -> _ tuple-> getData () [0]);
$ ey_y = floatval ($ ey -> _ tuple-> getData () [1]);
$ ey_z = floatval ($ ey -> _ tuple-> getData () [2]);
$ ey = nowy Math_Vector3 (tablica ($ ey_x, $ ey_y, $ ey_z));

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

// CALC D.
// zrób to wcześniej
$ d = floatval ($ d -> _ tuple-> getData () [0]);
// echo "d = $ d \ n";

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

# z wikipedii
#plug and chug przy użyciu powyższych wartości
$ 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);

# pokazano tylko jeden przypadek
$ z = sqrt (pow ($ DistA, 2) - pow ($ x, 2) - pow ($ y, 2));

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

#triPt to tablica z ECEF x, y, z punktu trilateracji
$ 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 = zmiennoprzecinkowy ($ triPt -> _ tuple-> getData () [0]);
$ triPt_y = zmiennoprzecinkowe ($ triPt -> _ tuple-> getData () [1]);
$ triPt_z = floatval ($ triPt -> _ tuple-> getData () [2]);


#convert powrót do lat / long z ECEF
#konwertuj na stopnie
$ lat = rad2deg (asin ($ triPt_z / $ earthR));
$ lon = rad2deg (atan2 ($ triPt_y, $ triPt_x));

echo $ lat. ','. $ lon;