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