Find k nearest neighbors on a sphere

Given a set $S$ of $N$ points on a sphere, and another point $P$ on the sphere, I want to find the $k$ points in $S$ that are the closest (Euclidean or great circle distance).

I'm willing to do a reasonable amount of pre-computation. The solution must be exact and efficient (faster than linear time).

• It looks like this question is easier to solve on the $\mathbb{R}^2$ plane. Have you thought about that, and do the answers translate easily to the sphere, and why (not)? For example, you can compute $A(v)=\{p | \forall_{w\in S} d(p,v) \leq d(p,w)\}$, which is a convex polygon defined by lines orthogonal to lines between $v$ and $w$, and use that together with space partitioning to find the nearest neighbour, and work from there. Oct 11 '15 at 21:37
• @LieuweVinkhuijzen the spherical topology matters, since the use case is a GIS problem. I have what I think will provide a working solution but am curious about prior art. Thanks! Oct 11 '15 at 22:25
• JohnJ, Lieuwe wasn't saying that the spherical topology is unimportant; instead, I think he is saying there are many standard techniques for solving this on $\mathbb{R}^2$, and I suspect he was recommending that you start by studying those standard approaches, and then look at whether they can be modified slightly to take into account that you're actually on a sphere rather than on $\mathbb{R}^2$.
– D.W.
Oct 12 '15 at 19:01

For instance, one approach is to use a $k$-d tree on on the surface of the sphere. You can express every point on the sphere using spherical coordinates: every point on the sphere has coordinate $(1,\theta,\phi)$. Thus, we have a 2-dimensional space with coordinates $(\theta,\phi)$. Now organize your points using a $k$-d tree, where here we are in $k=2$ dimensions. There are standard algorithms for nearest neighbor search in a $k$-d tree; heuristically, the expected running time is $O(\lg N)$.

You will have to make small modifications to the data structure to reflect that the coordinates "wrap around" modulo $2\pi$, but this is not hard. The key subroutine used in nearest neighbor search in a $k$-d tree is: given a point $P$ and a "rectangular" region $R$, find the distance from $P$ to the nearest point in $R$. In your case, the region $R$ is $[\theta_\ell,\theta_u] \times [\phi_\ell, \phi_u]$, i.e., the set of points $\{(1,\theta,\phi) : \theta_\ell \le \theta \le \theta_u, \phi_\ell \le \phi \le \phi_u\}$. It is easy to compute the distance from $P$ to the closest point in $R$. This will then let you use the standard algorithm for nearest neighbor search in a $k$-d tree.

Alternatively, instead of a $k$-d tree, you could use any other binary space partitioning tree, or you could look at metric trees, though I don't have any reason to expect them to be significantly better.

• I don't think this actually works, after all the k closest points in the unwrapped plane (even taking into account wrapping) are not necessarily the k closest points on the spherical surface, due to distortion. Nov 19 '21 at 14:59
• @matthias_buehlmann, OK. You may be right...
– D.W.
Nov 19 '21 at 21:45

Here are links to two different software packages that address your question. It may be worth studying each to see if the methods they employ satisfy your needs:

(1) Matlab GridSphere. "A geodesic grid is an even grid over the surface of a sphere. The algorithm is optimized for a grid generated by GridSphere and won't work on an arbitrary geodesic grid."

(2) DarkSkyApp sphere-knn .js. "provides fast nearest-neighbor lookups on a sphere....well-tested and works correctly regardless of where on the earth things are located."