Make a map of squares

Question

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[
   "http://code.google.com/apis/kml/documentation/us_states.kml",
   "Data"];
transform[s_] := StringTrim[s, Whitespace ~~ "(" ~~ ___ ~~ ")"]
polygons =
  Thread[transform[
     "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_] :=
 Apply[Plus,
   Flatten[First@poly Map[({1, -1} Reverse[#] &),
      RotateLeft[First@poly]]]]/2
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!

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_] :=

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

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

Timing[
 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]]) + 
     sl48/2}];
Graphics[Transpose[{colors, Rectangle @@@ rex2}]]