1
$\begingroup$

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?

$\endgroup$
4
$\begingroup$

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]]]}]

Mathematica graphics


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]}]

Mathematica graphics

$\endgroup$
  • $\begingroup$ 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? $\endgroup$ – ivanmp Mar 12 '15 at 0:45
  • $\begingroup$ 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"? $\endgroup$ – ivanmp Mar 12 '15 at 0:54
  • $\begingroup$ @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 $\endgroup$ – Dr. belisarius Mar 12 '15 at 1:34
  • $\begingroup$ Got it. I've generalized it now. Your solution is very nice, thank you! $\endgroup$ – ivanmp Mar 12 '15 at 11:38
1
$\begingroup$

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..)

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.