Wie gruppiere ich Nahpunkte mit GPS-Positionen?

Ich bin eine IT-Person, daher weiß ich nicht viel über Projektionen und so weiter. Ich hoffe, Sie können mir helfen.

Ich habe eine Anwendung für Android erstellt, die eine GPS-Position erfasst, sodass ich den Breiten- und Längengrad zu einem bestimmten Zeitpunkt habe. Ich möchte die Elemente, die nahe beieinander liegen, in Gruppen von Geländebereichen gleicher physischer Größe zusammen speichern. Das Problem ist, dass ich die Punkte vorher nicht kenne und sie von jeder Position auf der Welt kommen können .

Meine erste Idee (um das Problem etwas näher zu erläutern) war, die Dezimalstellen des Breiten- und Längengrads für die Gruppierung zu verwenden. Eine Gruppe sind beispielsweise die Positionen mit einem Breitengrad zwischen 35,123 und 35,124 und einem Längengrad zwischen 60,101 und 60,102. Wenn ich also eine Position wie lat = 35.1235647 und lon = 60.1012254598 erhalte, geht dieser Punkt an diese Gruppe.

Diese Lösung wäre für eine kartesische 2D-Darstellung in Ordnung, da ich Flächen von 0,001 Einheiten breit und hoch hätte; Da die Größe von 1 Längengrad in verschiedenen Breiten unterschiedlich ist, kann ich diesen Ansatz nicht verwenden.

Irgendeine Idee?

Ein schneller und schmutziger Weg verwendet eine rekursive sphärische Unterteilung . Beginnen Sie mit einer Triangulation der Erdoberfläche und teilen Sie jedes Dreieck rekursiv von einem Scheitelpunkt bis zur Mitte seiner längsten Seite. (Idealerweise teilen Sie das Dreieck in zwei Teile mit gleichem Durchmesser oder Teile mit gleicher Fläche, aber da diese eine umständliche Berechnung erfordern, teile ich die Seiten nur genau in zwei Hälften: Dies führt dazu, dass die verschiedenen Dreiecke letztendlich ein wenig unterschiedlich groß sind, aber das scheint für diese Anwendung nicht kritisch zu sein.)

Natürlich behalten Sie diese Unterteilung in einer Datenstruktur bei, die es ermöglicht, das Dreieck, in dem ein beliebiger Punkt liegt, schnell zu identifizieren. Ein Binärbaum (basierend auf den rekursiven Aufrufen) funktioniert gut: Jedes Mal, wenn ein Dreieck geteilt wird, wird der Baum am Knoten dieses Dreiecks geteilt. Daten bezüglich der Aufteilungsebene bleiben erhalten, sodass Sie schnell bestimmen können, auf welcher Seite der Ebene ein beliebiger Punkt liegt: Dies bestimmt, ob Sie den Baum nach links oder rechts hinunterfahren.

(Habe ich "Ebene" teilen gesagt? Ja - wenn Sie die Erdoberfläche als Kugel modellieren und geozentrische (x, y, z) Koordinaten verwenden, finden die meisten unserer Berechnungen in drei Dimensionen statt, wobei die Seiten der Dreiecke Teile von sind Schnittpunkte der Kugel mit Ebenen durch ihren Ursprung. Dies macht die Berechnungen schnell und einfach.)

Ich werde dies veranschaulichen, indem ich die Prozedur an einem Oktanten einer Kugel zeige. Die anderen sieben Oktanten werden auf die gleiche Weise verarbeitet. Ein solcher Oktant ist ein 90-90-90-Dreieck. In meinen Grafiken werde ich euklidische Dreiecke zeichnen, die sich über dieselben Ecken erstrecken: Sie sehen nicht sehr gut aus, bis sie klein werden, aber sie können einfach und schnell gezeichnet werden. Hier ist das euklidische Dreieck, das dem Oktanten entspricht: Es ist der Beginn der Prozedur.

Da alle Seiten gleich lang sind, wird eine zufällig als die "längste" ausgewählt und unterteilt:

Wiederholen Sie dies für jedes der neuen Dreiecke:

Nach n Schritten haben wir 2 ^ n Dreiecke. Hier ist die Situation nach 10 Schritten mit 1024 Dreiecken im Oktanten (und 8192 auf der Kugel insgesamt):

Zur weiteren Veranschaulichung habe ich einen zufälligen Punkt innerhalb dieses Oktanten generiert und den Unterteilungsbaum durchlaufen, bis die längste Seite des Dreiecks weniger als 0,05 Bogenmaß erreicht hat. Die (kartesischen) Dreiecke werden mit dem Sondenpunkt in Rot angezeigt.

Um die Position eines Punktes auf einen Breitengrad (ungefähr) zu beschränken, würden Sie übrigens feststellen, dass dies ungefähr 1/60 Radian ist und somit ungefähr (1/60) ^ 2 / (Pi / 2) = 1/6000 des Gesamtfläche. Da jede Unterteilung die Dreiecksgröße ungefähr halbiert, reichen etwa 13 bis 14 Unterteilungen des Oktanten aus. Das ist nicht viel Berechnung - wie wir weiter unten sehen werden -, was es effizient macht, den Baum überhaupt nicht zu speichern, sondern die Unterteilung im laufenden Betrieb durchzuführen. Zu Beginn würden Sie notieren, in welchem ​​Oktanten der Punkt liegt - das wird durch die Vorzeichen seiner drei Koordinaten bestimmt, die als dreistellige Binärzahl aufgezeichnet werden können - und bei jedem Schritt möchten Sie sich daran erinnern, ob der Punkt liegt links (0) oder rechts (1) des Dreiecks. Das ergibt eine weitere 14-stellige Binärzahl. Mit diesen Codes können Sie beliebige Punkte gruppieren.

(Wenn zwei Codes als tatsächliche Binärzahlen nahe beieinander liegen, sind die entsprechenden Punkte nahe beieinander. Punkte können jedoch immer noch nahe beieinander liegen und bemerkenswert unterschiedliche Codes aufweisen. Betrachten Sie beispielsweise zwei Punkte im Abstand von einem Meter vom Äquator: Ihre Codes müssen sich sogar unterscheiden vor dem Binärpunkt, weil sie sich in verschiedenen Oktanten befinden. So etwas ist bei jeder festen Aufteilung des Raums unvermeidlich.)

Ich habe Mathematica 8 verwendet , um dies zu implementieren: Sie können es so wie es ist oder als Pseudocode für eine Implementierung in Ihrer bevorzugten Programmierumgebung verwenden.

Bestimmen Sie, auf welcher Seite der Ebene 0-ab Punkt p liegt:

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

Verfeinern Sie das Dreieck abc basierend auf Punkt 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}] 
  ]

Die letzte Figur wurde gezeichnet, indem der Oktant angezeigt und darüber hinaus die folgende Liste als Satz von Polygonen gerendert wurde:

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]

NestWhileListrefineWendet wiederholt eine Operation ( ) an, während eine Bedingung vorliegt (das Dreieck ist groß) oder bis eine maximale Operationsanzahl erreicht wurde (16).

Um die vollständige Triangulation des Oktanten anzuzeigen, begann ich mit dem ersten Oktanten und wiederholte die Verfeinerung zehnmal. Dies beginnt mit einer geringfügigen Änderung von 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}}
  ]

Der Unterschied besteht darin , dass splitkehrt beide Hälften seines Eingangs Dreiecks statt das eine , bei der ein gegebener Punkt liegt. Die vollständige Triangulation wird erhalten, indem Folgendes wiederholt wird:

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

Zur Überprüfung habe ich ein Maß für die Größe jedes Dreiecks berechnet und mir den Bereich angesehen. (Diese "Größe" ist proportional zu der pyramidenförmigen Figur, die von jedem Dreieck und dem Mittelpunkt der Kugel begrenzt wird. Bei kleinen Dreiecken wie diesen ist diese Größe im Wesentlichen proportional zu ihrer Kugelfläche.)

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

{0,00523, 0,00739}

Daher variieren die Größen um etwa 25% gegenüber dem Durchschnitt. Dies erscheint sinnvoll, um einen annähernd einheitlichen Weg zur Gruppierung von Punkten zu erreichen.

Beim Scannen dieses Codes werden Sie keine Trigonometrie bemerken : Der einzige Ort, an dem er, wenn überhaupt, benötigt wird, ist das Hin- und Herkonvertieren zwischen sphärischen und kartesischen Koordinaten. Der Code projiziert auch nicht die Erdoberfläche auf eine Karte, wodurch die damit verbundenen Verzerrungen vermieden werden. Andernfalls werden nur die Mittelwertbildung ( Mean), der Satz von Pythagoras ( Norm) und eine 3 x 3-Determinante ( Det) verwendet, um die gesamte Arbeit zu erledigen. (Es gibt einige einfache Listenmanipulationsbefehle wie RotateLeftund Flattenauch zusammen mit einer Suche nach der längsten Seite jedes Dreiecks.)

Dies ist eine schwierige Frage, da Projektionen beim Übergang vom geografischen 3D-Koordinatensystem WGS84 zu einer flachen 2D-Projektion zu Verzerrungen führen. Auf globaler Ebene kann es irgendwo im System zu Verzerrungen kommen.

Ich denke, Ihre beste Wahl ist es, auf die Universal Transverse Mercator-Projektion zu projizieren . Soweit ich weiß, ist dies das Beste, was Sie einer globalen Projektion mit der geringsten Verzerrung erreichen können.

Wenn Sie die Anforderungen an die Definition von Gruppen in Bereichen mit genau derselben Größe sowie die Anforderungen an die Echtzeitverarbeitung lockern können, gibt es Clustering-Algorithmen wie DBSCAN und eine Familie von Derivaten, mit deren Hilfe Gruppierungen erstellt werden können.