Thanks for the many references to geodata libraries and implementations. But my question is more general: about the algorithms behind those implementations and their efficiency. To put the original question more abstractly: Given an arbitrary set of points on a sphere's surface, what is the most efficient algorithm to calculate a subset of the points which are within a certain distance over the sphere's surface to a certain (other) reference point on the sphere's surface? In addition a related problem: Given an arbitrary set of points on a sphere's surface, what is the most efficient algorithm to order this set according to the distance over the sphere's surface of each point to a certain (other) reference point on the sphere's surface?
The quite simple algorithm would be to calculate the distance of each point to the reference point and use these values for filtering or sorting respectively. But there should exist more sophisticated algorithms, like prefiltering or presorting the original set before applying more complex calculations or using other types of coordinate systems as latitude/longitude.
A trivial example would be if the filtering distance is more than half of the sphere circumference; in this case we do not need to calculate anything at all, because every point must be inside the distance. Or if we are using latitude/longitude notation and our reference point is on a pole, then we can just use the latitude for sorting without the need of distance calculation. One can easily think of some other special cases where specific alorithms may be used to speed up the operation.
Furthermore, I can think of a prefiltering algorithm for the latitude/longitude notation, where I calculate the coordinates of the four points directly north, south, west and east of the reference point in the given surface distance. Let us call them the "Noth point, "South point", "West point" and "East point". Then I exclude all points with coordinates west of the West point, east from the East point, north from the North point and south from the South point. Likewise I calculate the coordinates directly north-west, north-east, south-west, and south-east of the reference point in the given surface distance; all points within the minimum latitude/longitude range are within the given surface distance. Therefore I need only to perform the full calculation for the remaining points. (Special care must be taken if we cross the equator or the zero-meridian, of course.)
This is just a simple optimization that I came up with myself. I suppose, however, that a number of experts have already thought much more deeply about this group of problems. What interests me, is to learn about their solutions. An Internet search led me, for instance, to forum posts that recomment k-d trees to determine k nearest neighbors on a sphere. But I could not find a deeper analysis of this strategy regarding its computational efficiency anywhere.