Improving performace of `GeoDistance` and `NearestNeighborGraph`

Question

I am trying to construct a NearestNeighborGraph of the 3143 US counties, using their latitude and longitude positions.

counties = 
  EntityList[
   EntityClass["AdministrativeDivision", "USCountiesAllCounties"]];
countyPositions = EntityValue[counties, {"Latitude", "Longitude"}];

So far so good. However, computing NearestNeighborGraph[countyPositions, DistanceFunction -> GeoDistance] takes forever, or at least long enough that I aborted the computation because I got bored. It works as expected with a shorter list of counties, but is still pretty slow:

NearestNeighborGraph[Take[countyPositions, 10], 
  DistanceFunction -> GeoDistance] // AbsoluteTiming
(* {0.266751, graph omitted} *)

Also, if I just use latitude and longitude as Cartesian coordinates with the default EuclideanDistance it works just fine, and is considerably faster for the whole list of counties than using GeoDistance for just 10 of them!

NearestNeighborGraph[
  QuantityMagnitude /@ countyPositions] // AbsoluteTiming
(* {0.13108, graph omitted} *)

It makes sens that GeoDistance is slower than EuclideanDistance, but it's annoying that it's so much slower, and I would very much like to speed it up.

Answer 1

There are a few things that can be improved here:

1. Don't call {"Latitude", "Longitude"} separately in EntityValue. It's better to use the "Position" property, that returns a GeoPosition object. You also can use the entity class directly in EntityValue, without the need to send the 3143 county entities to WolframAlpha. Finally we wrap GeoPosition around the result to combine all positions into a single object only containing numbers. It also reduces the size of the object.

    countyPositions = GeoPosition[EntityValue[EntityClass["AdministrativeDivision", "USCountiesAllCounties"], "Position"]]

2. It is true that using GeoDistance is much slower than EuclideanDistance. The latter is a nearly trivial computation, while the former constructs a geodesic on an ellipsoid, essentially solving a differential equation to high precision per pair of locations. We can speed up computations by precomputing all distances in a single GeoDistance call:

    alldistances = GeoDistance[countyPositions, countyPositions]

That takes about 20 seconds for 10 million independent geodetic computations.

3. Extract all lonlat pairs from the GeoPosition above and setup numeric values for our own geodistance function:

    lonlats = countyPositions["LongitudeLatitude"];

    MapThread[Set, {Outer[geodistance, lonlats, lonlats, 1], QuantityMagnitude[alldistances, "Miles"]}, 2];

Now geodistance[lonlat1, lonlat2] is even faster than EuclideanDistance[lonlat1, lonlat2] because there is no computation involved. Just a downvalue lookup.

Now you can call

    NearestNeighborGraph[lonlats, DistanceFunction -> geodistance]

which will take about 40 seconds to produce its result. This is still slower than using EuclideanDistance, because NearestNeighborGraph is especially fast for a distance about which it knows all properties. Using an external distance function will be necessarily slower. Think of a geo computation involving points around the whole globe, with potential nearly-antipodal points, etc.