# How do I group near points with GPS positions?

I'm an IT person so I don't know too much about projections and so on, I hope you can help me.

I have made an application for Android that collects a GPS position, so I have the latitude and the longitude at a given time. I want to save together those elements which are near to each other, in groups of terrain areas of the same phisical size; the problem is that I don't know the points beforehand and they can come from any position in the world.

My first idea (to explain a little bit more the problem) was to use the decimals of the latitude and the longitude for grouping; for instance, one group will be the positions with latitude between 35.123 and 35.124, and the longitude between 60.101 and 60.102. So if I get a position like lat=35.1235647 and lon=60.1012254598, this point will go to that group.

This solution would be OK for a cartesian 2D representation, as I would have areas of 0.001 units wide and heigh; however, as the size of 1degree of longitude is different at different latitudes, I cannot use this approach.

Can you help me with this problem? Any idea?

-
Why can't you store the position as in, and then do the processing later? Also at the equator, 1 degree of longitude is approx 111 km, so 0.001 degree will be a little more than 1 Km. DO you really want your bins so Big? – Devdatta Tengshe Mar 11 '13 at 11:00
The 0.001 degree was just an example of my idea. Of course I will have to adapt it to the requirements. I cannot do a post processing as it is plannes to be a real time application, and I cannot say to my users "hold on until tomorrow, because I have to group the points" Thanks for the ideas anyway ;) – Kuu Mar 11 '13 at 14:51

One quick and dirty way uses a recursive spherical subdivision. Beginning with a triangulation of the earth's surface, recursively split each triangle from a vertex across to the middle of its longest side. (Ideally you will split the triangle into two equal-diameter parts or equal-area parts, but because those involve some fiddly calculation, I just split the sides exactly in half: this causes the various triangles eventually to vary a little in size, but that does not seem critical for this application.)

Of course you will maintain this subdivision in a data structure that allows quickly identifying the triangle in which any arbitrary point lies. A binary tree (based on the recursive calls) works nicely: each time a triangle is split, the tree is split at that triangle's node. Data concerning the splitting plane are retained, so that you can quickly determine on which side of the plane any arbitrary point lies: that will determine whether to travel left or right down the tree.

(Did I say splitting "plane"? Yes--if model the earth's surface as a sphere and use geocentric (x,y,z) coordinates, then most of our calculations take place in three dimensions, where sides of triangles are pieces of intersections of the sphere with planes through its origin. This makes the calculations fast and easy.)

I will illustrate by showing the procedure on one octant of a sphere; the other seven octants are processed in the same manner. Such an octant is a 90-90-90 triangle. In my graphics I will draw Euclidean triangles spanning the same corners: they don't look very good until they get small, but they can be easily and quickly drawn. Here is the Euclidean triangle corresponding to the octant: it is the start of the procedure.

As all sides have equal length, one is chosen randomly as the "longest" and subdivided:

Repeat this for each of the new triangles:

After n steps we will have 2^n triangles. Here is the situation after 10 steps, showing 1024 triangles in the octant (and 8192 on the sphere overall):

As a further illustration, I generated a random point within this octant and traveled the subdivision tree until the triangle's longest side reached less than 0.05 radians. The (cartesian) triangles are shown with the probe point in red.

Incidentally, to narrow a point's location to one degree of latitude (approximately), you would note that this is about 1/60 radian and so covers around (1/60)^2 / (Pi/2) = 1/6000 of the total surface. Since each subdivision approximately halves the triangle size, about 13 to 14 subdivisions of the octant will do the trick. That's not much computation--as we will see below--making it efficient not to store the tree at all, but to perform the subdivision on the fly. At the beginning you would note which octant the point lies in--that is determined by the signs of its three coordinates, which can be recorded as a three-digit binary number--and at each step you want to remember whether the point lies in the left (0) or right (1) of the triangle. That gives another 14-digit binary number. The full code for the illustrated point happens to be 111.10111110111100. You may use these codes to group arbitrary points.

(Generally, when two codes are close as actual binary numbers, the corresponding points are close; however, points can still be close and have remarkably different codes. Consider two points one meter apart separated the Equator, for instance: their codes must differ even before the binary point, because they are in different octants. This kind of thing is unavoidable with any fixed partitioning of space.)

I used Mathematica 8 to implement this: you may take it as-is or as pseudocode for an implementation in your favorite programming environment.

Determine what side of the plane 0-a-b point p lies on:

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

Refine triangle a-b-c based on 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}]
]
``````

The last figure was drawn by displaying the octant and, on top of that, by rendering the following list as a set of polygons:

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

`NestWhileList` repeatedly applies an operation (`refine`) while a condition pertains (the triangle is large) or until a maximum operation count has been reached (16).

To display the full triangulation of the octant, I began with the first octant and iterated the refinement ten times. This begins with a slight modification of `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}}
]
``````

The difference is that `split` returns both halves of its input triangle rather than the one in which a given point lies. The full triangulation is obtained by iterating this:

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

To check, I computed a measure of every triangle's size and looked at the range. (This "size" is proportional to the pyramid-shaped figure subtended by each triangle and the sphere's center; for small triangles like these, this size is essentially proportional to its spherical area.)

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

{0.00523, 0.00739}

Thus the sizes vary up or down by about 25% from their average: that seems reasonable for achieving an approximately uniform way to group points.

In scanning over this code, you will notice no trigonometry: the only place it will be needed, if at all, will be in converting back and forth between spherical and cartesian coordinates. The code also does not project the earth's surface onto any map, thereby avoiding the attendant distortions. Otherwise, it only uses averaging (`Mean`), the Pythagorean Theorem (`Norm`) and a 3 by 3 determinant (`Det`) to do all the work. (There are some simple list-manipulation commands like `RotateLeft` and `Flatten`, too, along with a search for the longest side of each triangle.)

-
+1 for a very thorough answer – Paul Hiemstra Mar 12 '13 at 13:33

This is a tough one, since projections introduce distortion when going from the 3D WGS84 geographic coordinate system to a flat 2D projection. At a global scale, you're bound to have distortions somewhere in the system.

I think your best bet is to project to the Universal Transverse Mercator projection. As far as I know this is the closest you can get to a global projection with the least amount of distortion.

If you can relax the requirements that groups have to be defined in areas of the exact same size as well as the requirement of real-time processing, there are clustering algorithms such as DBSCAN, and a family of derivatives, that can help produce groupings.

-
UTM is not a "global projection": the demonstration is that almost any pair of valid coordinates, such as (500000, 5000000), corresponds to at least 120 distinct widely-separated points in the UTM system. And, unfortunately, clustering algorithms do not meet the OP's needs of being able to group points in real time based solely on their location (and not on their proximity to other points). – whuber Mar 11 '13 at 22:20
@whuber re: "global projection" - right. That's why I said "the closest you can get to a global projection". If you know of a better projection system that's more appropriate please leave it in the comments and I'll edit my answer. Also, OP didn't have real-time requirement in the initial post. I'll edit my answer to take this into account. – Sean Barbeau Mar 12 '13 at 13:30
Sean, (1) My solution to the global projection problem is not to use one. There exists no global projection without singularities. (2) True, the real-time clarification appeared in a comment. The text following "my first idea" does do a good job, though, in emphasizing that this problem is one of partitioning the earth's surface rather than clustering a set of locations. That's the point I was (not very effectually) trying to get across in my earlier comment to you. – whuber Mar 12 '13 at 22:38