# Data structure to store sphere points (latitude,longitude) and retrieve all points within a distance

I have a set of thousands~millions of points on a sphere's surface, each with latitude, longitude.
I want to quickly get all points within a distance d of a particular sphere point latC, lonC.

• The points are mostly concentrated at places where humans go often.
• Some points have the exact same latitude and longitude.
• d is orders of magnitude smaller than the Earth's radius r, so using arc instead of distance is totally OK
• Getting a few extra points is OK (even a rectangle area instead of a circle area is better than nothing). Missing points is not OK
• No persistence needed.

What data structure would give these points quickly consistently, even if the set is large?

What I tried:

• Simple list: Checking all points is too slow.
• Quadtree on latitude, longitude: Does not find points around the poles.
• so you want to have a reference point and a distance and get all points within distance $d$ from this point as a list in one query lets say? i.e without filtering and searching points that match Commented Jan 8, 2016 at 9:59
• to get that list in $O(1)$ (constant) time, you would a data structure with $O(\max(d) n^2)$ memory which is quite large so some kind of filtering is needed, but such that is quite fast and optimal Commented Jan 8, 2016 at 10:01
• if the reference point was fixed then that would be easy since one would only have to sort the points by distance and just get the list of points up to a percentage $d/\max(d)$. But since the reference point is not fixed, this becomes dynamic (in a sense context-sensitive, by a slight abuse of terminology) Commented Jan 8, 2016 at 10:04
• Commented Jan 8, 2016 at 10:38
• Since you're only interested in small $d$, you could use multiple quadtrees to cover the entire map (transforming coordinates appropriately). Commented Jan 8, 2016 at 12:34

Geographic information systems use spherical triangular quadtrees to solve this problem. They are analogous to a quadtree, but defined over a sphere using triangular patches.

There are different variants with different bells and whistles, but you can refer to Quaternary Triangular Meshes (Dutton. "Encoding and Handling Geospatial Data with Hierarchical Triangular Meshes"), Sphere Quadtrees (Fekete. "Rendering and managing spherical data with sphere quadtree") and HTMs (Kunszt. "The Hierarchical Triangular Mesh") as approaches you can extrapolate from.

The basic idea is that you partition the sphere into coarse triangular patches, then subdivide those patches into smaller triangular patches recursively as needed. You can then perform basic queries like you would on a quadtree.

See Chen. "An Algorithm for the Generation of Voronoi Diagrams on the Sphere Based on QTM" for a survey of different techniques, and related problems.

Represent the points in cartesian coordinates (in 3D), and use square of distance (as it is cheaper to compute).

• In 3D! I see, that's doable then. I would still mean calculating for each point of the set, but each calculation would be faster indeed. Commented Jan 8, 2016 at 15:15
• @NicolasRaoul, you can e.g. sort the points by one of the coordinates to reduce it to a binary search. Commented Jan 8, 2016 at 16:02
• Just wanted to point out that this answer is simply incorrect. The constrained distance between two points on a sphere (i.e., their distance when you are forced to move along the surface of the sphere) is larger than their straight-line distance you would get from calculating their Euclidean distance. This is easier to see if you consider the simpler case of 2 points on a circle. Commented Jan 20, 2020 at 21:58
• @PMende, true. But the distance on the sphere is monotonic on linear distance. Commented Jan 29, 2020 at 17:07
• @PMende: Good remark, however the question mentions that d is orders of magnitude smaller than the Earth's radius r so it should be acceptable. Commented Feb 26, 2021 at 7:44

Use multiple quadtrees. Generate one quadtree based on lat-long. This has trouble at the poles (because the poles are "singularities", so pick a second lat-long-like coordinate system, but rotated 90 degrees, so the singularities are at two opposite points on the equator. Build a second quadtree based on that coordinate system. Now given any query, you can pick one of the two quadtrees and look it up in that quadtree: pick the quadtree that makes the query point as far away from the quadtree's singularities as possible.