...which raises the question: What is the most efficient way to store points on the sphere for lookup? Computationally? And in terms of storage?
1.) You can store lat/long, and use the Haversine formula, as you suggested. This requires trig functions, and has O(n) complexity; you need to iterate through all the points. You also have varying resolution over the surface, which makes bounding and early-outs a bit harder.
2.) A great many other coordinate charts also exist, and it's hard to say why you should choose one over the other without looking in detail at how the distance calculations are performed, etc.
3.) By using multiple charts -- e.g., a cube projection -- you can avoid issues with singularities, at the cost of branching. The complexity of distance calculations depends on the projection, but, without looking too carefully, my bet is that, in terms just of raw speed, cubemap vs. lat/long is probably a wash.
4.) Why use a coordinate chart at all, when you can use an embedding? If you store points in 3d, proximity calculations (since the points are on the sphere) just become a dot product. Much faster! It also opens up the possibility of, e.g. (if you will be doing many lookups but few insertions), storing indexes sorted along the three axes (or more!) to speed bounding-box (or more generally, sweep-'n-prune) calculations. Bins, bounding volume hierarchies, and the other standard tricks of computational geometry come into play. On the other hand, you're wasting a lot of codewords on points that don't actually lie on the sphere.
5.) Is there a more efficient use of codewords? Perhaps a (nearly-)-constant-resolution encoding scheme? If you start with the idea that a node in an octtree can be thought of as an octal number, you can see how you can encode points as real numbers in the interval [0, 1] -- e.g. "octal: .017135270661201") Of course, this still wastes codewords on points not on the sphere, so let's consider a refinement of this idea: At each level of the octree, discard any cube that does not intersect the sphere, and use arithmetic encoding, with the base varying between 8 and 2 depending on the number of cubes that intersect the sphere. This now seems like a (memory)-efficient way to encode points on the sphere -- but it is surely not computationally efficient. On the plus side, this same idea works for any manifold embedded in any Euclidean space, so at least it generalizes.
6.) Since #5 is a mapping from [0,1] to the sphere, one wonders if there are space-filling curves on the sphere. Of course there are -- e.g., the Hilbert curve in 2d, composed with any inverse coordinate chart. Not that this helps much!
I think my favorite of these is #5, but, practically, #1 or #4 are probably better choices.
So how do the real GiS systems do it?