Jak grupować w pobliżu punktów z pozycjami GPS?

Jestem informatykiem, więc nie wiem za dużo o prognozach i mam nadzieję, że możesz mi pomóc.

Zrobiłem aplikację na Androida, która zbiera pozycję GPS, więc mam szerokość i długość geograficzną w danym momencie. Chcę ocalić razem te elementy, które są blisko siebie, w grupach terenowych o tej samej wielkości fizycznej; Problem polega na tym , że nie znam punktów wcześniej i mogą pochodzić z dowolnej pozycji na świecie .

Moim pierwszym pomysłem (aby wyjaśnić nieco więcej problem) było użycie miejsc dziesiętnych szerokości i długości geograficznej do grupowania; na przykład jedną grupą będą pozycje o szerokości geograficznej między 35.123 a 35.124 i długości geograficznej między 60.101 a 60.102. Więc jeśli dostanę pozycję taką jak lat = 35,1235647 i lon = 60.1012254598, ten punkt przejdzie do tej grupy.

To rozwiązanie byłoby odpowiednie dla kartezjańskiej reprezentacji 2D, ponieważ miałbym obszary o szerokości 0,001 jednostki i wysokości; jednak ponieważ rozmiar 1 stopnia długości geograficznej jest różny na różnych szerokościach geograficznych, nie mogę zastosować tego podejścia.

Dowolny pomysł?

Jeden szybki i brudny sposób wykorzystuje rekurencyjny podział sferyczny . Zaczynając od triangulacji powierzchni ziemi, rekurencyjnie rozdzielaj każdy trójkąt od wierzchołka do środka jego najdłuższego boku. (Idealnie podzielisz trójkąt na dwie części o równej średnicy lub części o równej powierzchni, ale ponieważ wymagają one skomplikowanych obliczeń, po prostu podzielę boki dokładnie na pół: powoduje to, że różne trójkąty ostatecznie różnią się nieco wielkością, ale nie wydaje się to krytyczne dla tej aplikacji).

Oczywiście utrzymasz ten podział w strukturze danych, która pozwala szybko zidentyfikować trójkąt, w którym leży dowolny dowolny punkt. Drzewo binarne (oparte na wywołaniach rekurencyjnych) działa dobrze: za każdym razem, gdy trójkąt jest dzielony, drzewo jest dzielone w węźle tego trójkąta. Dane dotyczące płaszczyzny podziału są zachowywane, dzięki czemu można szybko określić, po której stronie płaszczyzny leży dowolny dowolny punkt: który określa, czy należy podróżować w lewo, czy w prawo po drzewie.

(Czy mówiłem o dzieleniu „płaszczyzny”? Tak - jeśli modelujemy powierzchnię ziemi jako kulę i używamy współrzędnych geocentrycznych (x, y, z), wówczas większość naszych obliczeń odbywa się w trzech wymiarach, gdzie boki trójkątów są kawałkami przecięcia sfery z płaszczyznami poprzez jej pochodzenie. Dzięki temu obliczenia są szybkie i łatwe).

Zilustruję to pokazując procedurę na jednej oktanie kuli; pozostałe siedem oktantów jest przetwarzanych w ten sam sposób. Taki oktan to trójkąt 90–90–90. W mojej grafice narysuję trójkąty euklidesowe obejmujące te same rogi: nie wyglądają zbyt dobrze, dopóki nie staną się małe, ale można je łatwo i szybko narysować. Oto trójkąt euklidesowy odpowiadający oktanowi: jest to początek procedury.

Ponieważ wszystkie boki mają równą długość, jedna jest wybierana losowo jako „najdłuższa” i dzieli się na:

Powtórz to dla każdego z nowych trójkątów:

Po n krokach będziemy mieli 2 ^ n trójkątów. Oto sytuacja po 10 krokach, pokazująca 1024 trójkąty w oktanie (i 8192 na kuli ogółem):

Jako kolejną ilustrację wygenerowałem losowy punkt w tej oktanie i przemierzyłem drzewo podziału, aż najdłuższy bok trójkąta osiągnie mniej niż 0,05 radianów. Trójkąty (kartezjańskie) są oznaczone kolorem punktu próbkowania na czerwono.

Nawiasem mówiąc, aby zawęzić położenie punktu do jednego stopnia szerokości geograficznej (w przybliżeniu), należy zauważyć, że jest to około 1/60 radianów, a więc obejmuje około (1/60) ^ 2 / (Pi / 2) = 1/6000 z powierzchnia całkowita. Ponieważ każdy poddział dzieli w przybliżeniu o połowę rozmiar trójkąta, około 13 do 14 pododdziałów oktanu załatwi sprawę. To niewiele obliczeń - jak zobaczymy poniżej - czyniąc efektywnym nie przechowywanie drzewa w ogóle, ale dokonywanie podziału w locie. Na początku należy zauważyć, w której ósemce znajduje się punkt - który jest określony przez znaki jego trzech współrzędnych, które można zapisać jako trzycyfrową liczbę dwójkową - i na każdym kroku chcesz pamiętać, czy punkt leży w lewej (0) lub prawej (1) trójkąta. To daje kolejną 14-cyfrową liczbę binarną. Możesz użyć tych kodów do grupowania dowolnych punktów.

(Zasadniczo, gdy dwa kody są bliskie jak rzeczywiste liczby binarne, odpowiednie punkty są bliskie; jednak punkty mogą nadal być bliskie i mieć wyjątkowo różne kody. Weźmy na przykład dwa punkty oddzielone równikiem o jeden metr: ich kody muszą się różnić nawet przed punktem binarnym, ponieważ są w różnych oktantach. Tego rodzaju rzeczy nie da się uniknąć przy żadnym stałym podziale przestrzeni.)

Użyłem Mathematica 8, aby to zaimplementować: możesz wziąć to tak, jak jest lub jako pseudokod do implementacji w twoim ulubionym środowisku programistycznym.

Określ, na której stronie płaszczyzny punkt 0-ab p leży:

side[p_, {a_, b_}] := If[Det[{p, a, b}] >=  0, left, right];

Udoskonal trójkąt abc na podstawie punktu p.

refine[p_, {a_, b_, c_}] := Block[{sides, x, y, z, m},
  sides = Norm /@ {b - c, c - a, a - b} // N;
  {x, y, z} = RotateLeft[{a, b, c}, First[Position[sides, Max[sides]]] - 1];
  m = Normalize[Mean[{y, z}]];
  If[side[p, {x, m}] === right, {y, m, x}, {x, m, z}] 
  ]

Ostatnia cyfra została narysowana przez wyświetlenie oktanta, a ponadto przez renderowanie poniższej listy jako zestawu wielokątów:

p = Normalize@RandomReal[NormalDistribution[0, 1], 3]        (* Random point *)
{a, b, c} = IdentityMatrix[3] . DiagonalMatrix[Sign[p]] // N (* First octant *)
NestWhileList[refine[p, #] &, {a, b, c}, Norm[#[[1]] - #[[2]]] >= 0.05 &, 1, 16]

NestWhileListwielokrotnie stosuje operację ( refine), dopóki spełniony jest warunek (trójkąt jest duży) lub do momentu osiągnięcia maksymalnej liczby operacji (16).

Aby wyświetlić pełną triangulację oktanu, zacząłem od pierwszego oktanu i powtórzyłem wyrafinowanie dziesięć razy. Zaczyna się to od niewielkiej modyfikacji refine:

split[{a_, b_, c_}] := Module[{sides, x, y, z, m},
  sides = Norm /@ {b - c, c - a, a - b} // N;
  {x, y, z} = RotateLeft[{a, b, c}, First[Position[sides, Max[sides]]] - 1];
  m = Normalize[Mean[{y, z}]];
  {{y, m, x}, {x, m, z}}
  ]

Różnica polega na tym, że splitzwraca obie połówki trójkąta wejściowego, a nie tę, w której leży dany punkt. Pełna triangulacja jest uzyskiwana przez iterację:

triangles = NestList[Flatten[split /@ #, 1] &, {IdentityMatrix[3] // N}, 10];

Aby to sprawdzić, obliczyłem miarę każdego trójkąta i spojrzałem na zakres. (Ten „rozmiar” jest proporcjonalny do figury w kształcie piramidy, znajdującej się w każdym trójkącie i środku kuli; w przypadku takich małych trójkątów rozmiar ten jest zasadniczo proporcjonalny do jego sferycznego obszaru).

Through[{Min, Max}[Map[Round[Det[#], 0.00001] &, triangles[[10]] // N, {1}]]]

{0,00523, 0,00739}

Zatem rozmiary różnią się w górę lub w dół o około 25% od ich średniej: wydaje się to rozsądne, aby osiągnąć w przybliżeniu jednolity sposób grupowania punktów.

Podczas skanowania tego kodu nie zauważysz żadnej trygonometrii : jedynym miejscem, w którym będzie ono potrzebne, będzie w ogóle konwersja w obie strony między współrzędnymi sferycznymi i kartezjańskimi. Kod nie rzutuje również powierzchni ziemi na żadną mapę, co pozwala uniknąć towarzyszących jej zniekształceń. W przeciwnym razie do wykonania całej pracy używa tylko uśredniania ( Mean), twierdzenia Pitagorasa ( Norm) i wyznacznika 3 na 3 ( Det). (Istnieje kilka prostych poleceń list-manipulacja, jak RotateLefti Flattenteż, wraz z poszukiwaniem najdłuższego boku każdego trójkąta).

Jest to trudne, ponieważ projekcje wprowadzają zniekształcenie podczas przechodzenia z układu współrzędnych geograficznych 3D WGS84 do płaskiej projekcji 2D. W skali globalnej na pewno gdzieś w systemie występują zniekształcenia.

Myślę, że najlepszym rozwiązaniem jest rzutowanie na projekcję Universal Transverse Mercator . O ile wiem, jest to najbliżej globalnej projekcji z najmniejszym zniekształceniem.

Jeśli możesz złagodzić wymagania, które muszą być zdefiniowane w obszarach o dokładnie takiej samej wielkości, a także wymagania dotyczące przetwarzania w czasie rzeczywistym, istnieją algorytmy klastrowania, takie jak DBSCAN , i rodzina pochodnych, które mogą pomóc w tworzeniu grup.