Comment regrouper des points proches avec des positions GPS?

Je suis un informaticien donc je ne connais pas trop les projections et ainsi de suite, j'espère que vous pourrez m'aider.

J'ai fait une application pour Android qui recueille une position GPS, j'ai donc la latitude et la longitude à un moment donné. Je veux sauvegarder ensemble les éléments qui sont proches les uns des autres, dans des groupes de zones de terrain de même taille physique; le problème est que je ne connais pas les points à l'avance et qu'ils peuvent provenir de n'importe quelle position dans le monde .

Ma première idée (pour expliquer un peu plus le problème) était d'utiliser les décimales de la latitude et de la longitude pour le regroupement; par exemple, un groupe sera les positions avec une latitude comprise entre 35,123 et 35,124 et la longitude entre 60,101 et 60,102. Donc, si j'obtiens une position comme lat = 35.1235647 et lon = 60.1012254598, ce point ira à ce groupe.

Cette solution serait OK pour une représentation 2D cartésienne, car j'aurais des zones de 0,001 unités de large et de haut; cependant, comme la taille de 1 degré de longitude est différente à différentes latitudes, je ne peux pas utiliser cette approche.

Une idée?

Une manière rapide et sale utilise une subdivision sphérique récursive . Commençant par une triangulation de la surface de la terre, divisez récursivement chaque triangle d'un sommet à travers au milieu de son côté le plus long. (Idéalement, vous diviserez le triangle en deux parties de diamètre égal ou des parties de surface égale, mais parce que celles-ci impliquent des calculs délicats, je divise juste les côtés exactement en deux: cela fait que les différents triangles varient finalement un peu en taille, mais cela ne semble pas critique pour cette application.)

Bien sûr, vous conserverez cette subdivision dans une structure de données qui permet d'identifier rapidement le triangle dans lequel se situe tout point arbitraire. Un arbre binaire (basé sur les appels récursifs) fonctionne bien: chaque fois qu'un triangle est divisé, l'arbre est divisé au nœud de ce triangle. Les données concernant le plan de fractionnement sont conservées, afin que vous puissiez déterminer rapidement de quel côté du plan se situe un point quelconque: cela déterminera si vous devez vous déplacer à gauche ou à droite dans l'arbre.

(Ai-je dit diviser "plan"? Oui - si vous modélisez la surface de la terre comme une sphère et utilisez des coordonnées géocentriques (x, y, z), alors la plupart de nos calculs ont lieu en trois dimensions, où les côtés des triangles sont des morceaux de intersections de la sphère avec des plans passant par son origine. Cela rend les calculs rapides et faciles.)

Je vais illustrer en montrant la procédure sur un octant d'une sphère; les sept autres octants sont traités de la même manière. Un tel octant est un triangle 90-90-90. Dans mes graphiques, je vais dessiner des triangles euclidiens couvrant les mêmes coins: ils ne semblent pas très bons jusqu'à ce qu'ils deviennent petits, mais ils peuvent être facilement et rapidement dessinés. Voici le triangle euclidien correspondant à l'octant: c'est le début de la procédure.

Comme tous les côtés ont une longueur égale, l'un est choisi au hasard comme le "plus long" et subdivisé:

Répétez cette opération pour chacun des nouveaux triangles:

Après n étapes, nous aurons 2 ^ n triangles. Voici la situation après 10 étapes, montrant 1024 triangles dans l'octant (et 8192 sur la sphère globale):

Comme illustration supplémentaire, j'ai généré un point aléatoire dans cet octant et parcouru l'arbre de subdivision jusqu'à ce que le côté le plus long du triangle atteigne moins de 0,05 radians. Les triangles (cartésiens) sont représentés avec le point de sonde en rouge.

Soit dit en passant, pour réduire la position d'un point à un degré de latitude (environ), vous remarquerez qu'il s'agit d'environ 1/60 radian et couvre donc environ (1/60) ^ 2 / (Pi / 2) = 1/6000 de la surface totale. Étant donné que chaque subdivision divise environ par deux la taille du triangle, environ 13 à 14 subdivisions de l'octant feront l'affaire. Ce n'est pas beaucoup de calcul - comme nous le verrons ci-dessous - ce qui rend efficace de ne pas stocker l'arbre du tout, mais d'effectuer la subdivision à la volée. Au début, vous noterez dans quel octant se situe le point - déterminé par les signes de ses trois coordonnées, qui peuvent être enregistrées sous la forme d'un nombre binaire à trois chiffres - et à chaque étape, vous voulez vous souvenir si le point se trouve à gauche (0) ou à droite (1) du triangle. Cela donne un autre nombre binaire à 14 chiffres. Vous pouvez utiliser ces codes pour regrouper des points arbitraires.

(Généralement, lorsque deux codes sont proches en tant que nombres binaires réels, les points correspondants sont proches; cependant, les points peuvent toujours être proches et avoir des codes remarquablement différents. Considérez deux points à un mètre de distance séparés de l'équateur, par exemple: leurs codes doivent différer même avant le point binaire, car ils sont dans des octants différents. Ce genre de chose est inévitable avec tout partitionnement fixe de l'espace.)

J'ai utilisé Mathematica 8 pour implémenter ceci: vous pouvez le prendre tel quel ou comme pseudocode pour une implémentation dans votre environnement de programmation préféré.

Déterminez de quel côté du point 0-ab plan p se trouve:

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

Affiner le triangle abc en fonction du point 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}] 
  ]

La dernière figure a été dessinée en affichant l'octant et, en plus de cela, en rendant la liste suivante comme un ensemble de polygones:

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]

NestWhileListapplique de façon répétée une opération ( refine) tant qu'une condition se produit (le triangle est grand) ou jusqu'à ce qu'un nombre maximal d'opérations soit atteint (16).

Pour afficher la triangulation complète de l'octant, j'ai commencé par le premier octant et j'ai répété le raffinement dix fois. Cela commence par une légère modification de 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}}
  ]

La différence est que splitrenvoie les deux moitiés de son triangle d'entrée plutôt que celle dans laquelle se trouve un point donné. La triangulation complète est obtenue en itérant ceci:

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

Pour vérifier, j'ai calculé une mesure de la taille de chaque triangle et regardé la plage. (Cette "taille" est proportionnelle à la figure en forme de pyramide sous-tendue par chaque triangle et le centre de la sphère; pour de petits triangles comme ceux-ci, cette taille est essentiellement proportionnelle à sa zone sphérique.)

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

{0,00523, 0,00739}

Ainsi, les tailles varient d'environ 25% vers le haut ou vers le bas par rapport à leur moyenne: cela semble raisonnable pour obtenir une manière approximativement uniforme de grouper les points.

En parcourant ce code, vous ne remarquerez aucune trigonométrie : le seul endroit où il sera nécessaire, voire pas du tout, sera la conversion dans les deux sens entre les coordonnées sphériques et cartésiennes. Le code ne projette pas non plus la surface de la Terre sur aucune carte, évitant ainsi les distorsions qui en découlent. Sinon, il utilise uniquement la moyenne ( Mean), le théorème de Pythagore ( Norm) et un déterminant 3 par 3 ( Det) pour faire tout le travail. (Il existe quelques commandes de manipulation de liste simples comme RotateLeftet Flattenaussi avec une recherche du côté le plus long de chaque triangle.)

Ceci est difficile, car les projections introduisent une distorsion lors du passage du système de coordonnées géographiques 3D WGS84 à une projection 2D plate. À l'échelle mondiale, vous êtes lié à des distorsions quelque part dans le système.

Je pense que votre meilleur pari est de projeter sur la projection Universal Transverse Mercator . Pour autant que je sache, c'est le plus proche d'une projection globale avec le moins de distorsion.

Si vous pouvez assouplir les exigences selon lesquelles les groupes doivent être définis dans des zones de la même taille ainsi que l'exigence de traitement en temps réel, il existe des algorithmes de clustering tels que DBSCAN et une famille de dérivés, qui peuvent aider à produire des regroupements.