Matrix of boundary distance between countries

I’m trying to generate a matrix of the boundary distances between all countries.

First thing I tried was GeoDistance:

table = ParallelTable[
QuantityMagnitude@
GeoDistance[x, y, DistanceFunction -> “Boundary”,
UnitSystem -> “Metric”],
{x, countries},
{y, countries}
];

But this is way too slow. I tried preloading the CountryData paclet, but it didn’t make any difference. GeoDistance is calling Wolfram’s servers for each and every pair of countries.

CountryData[All, “Preload”]
RebuildPacletData[]

I also tried using GeoPosition with CountryData[COUNTRY, “Coordinates”], but I can’t figure out how to make this work with GeoDistance for “Boundary” distance.

Additionally I’ve looked at geographic questions here, but I haven’t found any that is similar enough to mine to give me directions on this.

Is there a way I can do this computation offline and fast?

=================

=================

2 Answers
2

=================

The following is a sketch, where I considered only two countries and only one boundary curve description for each of them. Generalization for multiple boundaries at the end.

slv = CountryData[“El Salvador”, “Coordinates”][[1]];
arg = CountryData[“Argentina”, “Coordinates”][[2]];
slvInt = Interpolation /@ Transpose@slv
argInt = Interpolation /@ Transpose@arg
sol=Timing@Quiet@
NMinimize[{GeoDistance @@ Through /@ {argInt[r], slvInt[s]},
1 <= s <= 161 && 1 <= r <= 955}, {r, s}, Method -> #] & /@ {Automatic, “SimulatedAnnealing”,
“DifferentialEvolution”, “NelderMead”, “RandomSearch”}

Timings and results on a very slow machine:

{{53.609375, {4.53166*10^6, {r -> 930.,s -> 106.882}}},
{12.359375, {4.56584*10^6, {r -> 930.472,s -> 120.92}}},
{54.937500, {4.53166*10^6, {r -> 930., s -> 106.882}}},
{12.218750, {4.53401*10^6, {r -> 930., s -> 106.159}}},
{16.968750, {4.56692*10^6, {r -> 930., s -> 128.242}}}}

Testing:

Graphics[{Gray, EdgeForm[Black],
CountryData[#, “Polygon”] & /@ CountryData[“SouthAmerica”],
CountryData[#, “Polygon”] & /@ CountryData[“NorthAmerica”], Thick,
Red, Line[ Reverse /@ Through /@ {argInt[r], slvInt[s]} /. sol[[1, 2, 2]]]}]

Generalization for multiple boundaries

Needs[“DifferentialEquations`InterpolatingFunctionAnatomy`”];
ifd[x_] := InterpolatingFunctionDomain[x][[1, 2]];
gtb = CountryData[“GreatBritain”, “Coordinates”];
ice = CountryData[“Iceland”, “Coordinates”];
borders = Map[Interpolation /@ (Transpose@#) &, #] & /@ {gtb, ice};
tb = Tuples[borders];

sol = (First@
Sort[MapIndexed[{Quiet@
NMinimize[{GeoDistance @@ Through /@ {#[[1]][r], #[[2]][s]},
1 <= s <= ifd[#[[2, 1]]] && 1 <= r <= ifd[#[[1, 1]]]}, {r, s}, Method -> “SimulatedAnnealing”] &@#1, #2} &, tb]])

{{distance, rule}, {idx}} = sol

Graphics[{Gray, EdgeForm[Black],
CountryData[#, “Polygon”] & /@ CountryData[“Europe”], Thick, Red,
Line[Reverse /@ Through /@ {tb[[idx,1]][r], tb[[idx,2]][s]} /. rule]}]

  

 

Thank you for your help! This is actually along the lines of what I wanted to do (but didn’t know how). Unfortunately, I can’t reproduce the results: bit.ly/1MsRnDE. Any idea what happened?
– ivanmp
Mar 12 ’15 at 0:45

  

 

Ignore last comment. Turns out all I had to do was to add QuantityMagnitude to the NMinimize to turn the results from miles into just numbers. One last thing: What did you mean by “only one boundary curve description for each of them”?
– ivanmp
Mar 12 ’15 at 0:54

  

 

@ivanmp Each separate territory in a country (for example an island) gives you a polygon in CountryData[“xxx”, “Coordinates”]. In my example I considered only one polygon for each country to abbreviate the list mapping coding and make the idea easier to follow
– Dr. belisarius
Mar 12 ’15 at 1:34

  

 

Got it. I’ve generalized it now. Your solution is very nice, thank you!
– ivanmp
Mar 12 ’15 at 11:38

fast approximation..

slv = Flatten[CountryData[“El Salvador”, “Coordinates”],1];
arg = Flatten[CountryData[“Argentina”, “Coordinates”],1];
(f = Nearest[slv];
Min[GeoDistance[First@f@# , #] & /@ arg]) // Timing

{2.012413, 4.53175*10^6}

This finding the minimum Norm of the angles.. so will fail in cases where the min distance is across the date line, etc.

( You can supply GeoDistance as a DistanceFunction to Nearest but its very slow..)