I'm trying to use Mathematica to make a map similar to the one here


i.e. a US map where each state is represented by a square of a given size, located as close as possible to its true geographical position while not overlapping with other squares. I get the data for the squares as follows:

usa = Import[
transform[s_] := StringTrim[s, Whitespace ~~ "(" ~~ ___ ~~ ")"]
polygons =
     "PlacemarkNames" /. usa[[1]]] -> ("Geometry" /. usa[[1]])];
(* Remove Alaska and Hawai *)
polygons =
  Join[{polygons[[1]]}, polygons[[3 ;; 10]], polygons[[12 ;;]]];
centers = First /@ polygons[[All, 2, 1]];
area[poly_] :=
   Flatten[First@poly Map[({1, -1} Reverse[#] &),
areas = Table[Total[area /@ polygons[[i, 2, 2 ;;]]], {i, 1, 48}];
sizes = Sqrt /@ areas;

which uses the true area of the state as square area, which is enough for now. I can plot them:

Graphics[Table[{RGBColor[RandomReal[], RandomReal[], RandomReal[]], 
   Rectangle[{centers[[i, 1]] - sizes[[i]]/2, 
     centers[[i, 2]] - sizes[[i]]/2},
    {centers[[i, 1]] + sizes[[i]]/2, 
     centers[[i, 2]] + sizes[[i]]/2}]}, {i, 1, 48}]]


Starting now, I demonstrate how far I've managed to go on my own. I want to move the squares to non-overlapping positions. I try to minimize the following score function:

dist[c1_, c2_, s1_, s2_] := Module[{dx, dy},
   dx = Abs[c1[1] - c2[1]]; dy = Abs[c1[2] - c2[2]];
   Max[0, dx - s1 - s2] + Max[0, dy - s1 - s2]
res = FindMinimum[
  Sum[Sum[dist[c[i], c[j], sizes[[i]], sizes[[j]]], {j, i + 1, 48}], {i, 1, 48}], 
  Flatten[Table[{c[i][j], centers[[i, j]]}, {i, 1, 48}, {j, 1, 2}], 1]]

but plotting this gives terrible results:

Graphics[Table[{RGBColor[RandomReal[], RandomReal[], RandomReal[]], 
    Rectangle[{c[i][1] - sizes[[i]]/2, c[i][2] - sizes[[i]]/2},
     {c[i][1] + sizes[[i]]/2, c[i][2] + sizes[[i]]/2}]}, {i, 1, 48}] /. res[[2]]]

I have not been able to see what is the issue with my code, because the score function clearly shouldn't be zero for this particular solution, but I am at loss to understand why. Any hint about that, or advice on tackling the problem in a different way, is welcome!

Doesn't work

  • 4
    $\begingroup$ Looks like a job for belisarius' answer to spread out rectangles to make them non-overlapping. Perhaps also subject it to a shape constraint (like in Heike's answer to the word cloud question) $\endgroup$
    – rm -rf
    Commented Sep 10, 2012 at 21:00
  • $\begingroup$ @R.M yeah, I saw the word-cloud question (and answers). I'd prefer something that clearly uses Mathematica's optimizers rather than a specially crafted algorithm, as it is also intended as a demonstration of Mathematica's power to solve real-life tasks… $\endgroup$
    – F'x
    Commented Sep 10, 2012 at 21:03
  • $\begingroup$ Did you see the link to belisarius' answer? There's a link to the code in the comments. Looks like it's a pretty straightforward application of that... I'll leave this for him to answer it for after all, he did all the hard work :) $\endgroup$
    – rm -rf
    Commented Sep 10, 2012 at 21:06
  • $\begingroup$ I had the same buttonbar problem tonight and I saw this complaint elsewhere too. $\endgroup$ Commented Sep 10, 2012 at 21:09
  • 2
    $\begingroup$ @R.M The code in that answer is crap, because it was done (almost) procedurally on purpose. I'll see if I can remember what I did and recast it functionally $\endgroup$ Commented Sep 11, 2012 at 5:35

1 Answer 1


First set up the lower 48 lists, along with a set of weighted offset vectors. I'll use these for initial points in a constrained optimization. They help insofar as they can be adjusted to satisfy, or almost satisfy, the inequality constraints (which, inexplicably, I call "penalty" below).

sl48 = Take[sizes, 48];
cenl48 = Take[centers, 48];
middle = Total[cenl48]/48;
wmean = sl48.cenl48/Total[sl48]
offsets = Map[# - wmean &, cenl48];

Our function will measure how far the new centers are from the initial ones. This is needed because otherwise we tend to get a result that is a tight fit of weighted squares, but bears little resemblance to the USA map.

f[centers : {{_?NumberQ, _?NumberQ} ..}, 
  orig : {{_?NumberQ, _?NumberQ} ..}] := With[
  {diffs = centers - orig},
  Total@Map[#.# &, diffs]]

overlap[c1_, c2_, s1_, s2_] := Min[(s1 + s2)/2 - Abs[c1 - c2]] <= 0

penalty[centers_, wts_] :=

   overlap[centers[[i]], centers[[j]], wts[[i]], wts[[j]]], {i, 
    47}, {j, i + 1, 48}]]

cvars = Array[c, {48, 2}];

 res = FindMinimum[{f[cvars, cenl48], penalty[cvars, sizes]}, 
    Transpose[{Flatten[cvars], Flatten[cenl48 + offsets]}], 
    MaxIterations -> 50, Method -> "InteriorPoint"];]

(* Out[344]= {136.37, Null} *)

We can plot this as follows.

colors = Table[
   RGBColor[RandomReal[], RandomReal[], RandomReal[]], {48}];
rex2 = Transpose[{(cvars /. res[[2]]) - sl48/2, (cvars /. res[[2]]) + 
Graphics[Transpose[{colors, Rectangle @@@ rex2}]]

enter image description here

  • 2
    $\begingroup$ I guess I should mention that this can probably be done as well or better, and much faster, if recast as a linear programming problem using Manhattan distances rather than Euclidean in the objective function (and slightly changing the constraints an explicitly linear formulation). $\endgroup$ Commented Sep 12, 2012 at 15:20

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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