Can this geo-3D histogram of USA population be improved?

Questions
,
1

I am looking for an optimal code to build a 3D histogram of USA population with bars bounded by state borders. Here is my initial version. Can anyone suggest improvements?

(* get the states *)
divisions = 
  EntityValue[Entity["AdministrativeDivision", {_, "UnitedStates"}], 
   "Entities"];

(* get polygons of borders *)
dat = EntityValue[
    divisions, {"Population", "Polygon"}] /. {GeoPosition -> Identity,
     Quantity[x_, _] -> x};

(* some arbitrary rescaling to improve relative height perception *)
pop = Rescale[(# - Min[#]) &@Log[dat[[All, 1]]] // N];

(* plot constants of population of regions of polygons *)
polygs = Plot3D[#1, {x, y} \[Element] #2, Mesh -> None, Filling -> 0, 
     ColorFunction -> "Rainbow", ColorFunctionScaling -> False] & @@@ 
   Transpose[{pop, dat[[All, 2]]}];

(* combine all *)
Show[polygs, PlotRange -> {{23, 50}, {-60, -130}, All}, 
 BoxRatios -> {27, 70, 50}, ImageSize -> 800, Boxed -> False, 
 Axes -> False]

enter image description here

enter image description here

2

Dear Vitaliy,

very nice indeed! I wanted to have a look at that but interestingly my system:

“10.0 for Mac OS X x86 (64-bit) (June 29, 2014)”

reproducibly crashes during the execution of the last command. I can send you the crash report.

Best wishes,
Marco

3

Thanks for the kind words and sorry for the crash. Please do send the report. I wonder if anyone else experience this.

4

Dear Vitaliy,

I have submitted the report. Also I tried four Macs with slightly different specifications, all on Mavericks, and they all crash at the same command.

Cheers,
Marco

5

Duplicated the crash on a MacOS 10.9.4 machine. Ran fine on MacOS 10.8.5 and Windows 7.

Mixed results on Linux.

Problem report filed.

6

Hello,
Vitaliy’s code doesn’t crash my MacBook MacOS 10.9.4 (16GB Graphics NVIDIA GeForce GT 750M 2048 MB)

I tried a different angle on Vitaliy’s code, but I don’t think it is an improvement. I tried a 2d version of the data representation because some of the states were getting hidden behind the others.

This is the same as Vitaliy’s code except I am grabbing the position and scaling a bit differently:

divisions = 
  EntityValue[Entity["AdministrativeDivision", {_, "UnitedStates"}], 
   "Entities"];
dat = EntityValue[
    divisions, {"Population", "Polygon", 
     "Position"}] /. {GeoPosition -> Identity, Quantity[x_, _] -> x};
dat[[All, 1]] = 
 Rescale[dat[[All, 1]], Through[{Min, Max}[dat[[All, 1]]]], {0.1, 1}]

Here is the 2D graphic where I have scaled the polygon of each state with its population and used a color scaling as well:

sizeScale[{population_, polygon_, center_}] := 
 With[{polygons =  polygon[[1]]},(*Print[polygons];*)
  polygons /. {x_?NumericQ, y_?NumericQ} :> 
    Reverse@(center + population ({x, y} - center))]

Show[
 Graphics[
  {FaceForm[], EdgeForm[Black], 
   dat[[All, 2]] /. p : {x_?NumericQ, y_?NumericQ} :> Reverse[p] }],
 
 Graphics[
  Transpose[{ColorData["DarkRainbow"] /@ dat[[All, 1]], 
    Polygon /@ (sizeScale /@ dat)}]], ImageSize -> Large
 ]

states.png

As you can see, there is a strange outlier near the international date line. It looks like the Philippines, but that isn’t one of the Entities in divisions (???)

7

I like Craigs idea of visualization very much. A lot of white space is generated by presence of islands. Here is the central pieces to enjoy the approach fully:

2014-08-24_14-14-30.png
2014-08-24_14-14-30.png1107×489 96.5 KB

8

Here’s another approach (it’s bit complicated but faster..) :

(* divide polygon pts to clean up artificials when polygon has holes *)
FindContourBreaks[pts_List] := 
  Module[{i, lines, breaks = {}}, 
   lines = {pts[[#[[1]]]], pts[[#[[2]]]]} & /@ 
     Partition[RotateLeft[Flatten[{#, #} & /@ Range[Length[pts]], 1]],
       2];
   Position[lines, 
     Alternatives @@ 
      Intersection[{lines[[All, 2]], lines[[All, 1]]} // Transpose, 
       lines]] // Flatten
   ];

FindContourBreak[pts_List] := 
  Module[{breaks, ranges}, breaks = FindContourBreaks[pts];
   ranges = 
    Partition[
     RotateLeft[Join[{1, 1}, Flatten[{#, # + 1} & /@ breaks]]], 2];
   ranges = Drop[ranges, -1];
   DeleteCases[Range @@@ ranges, x_ /; Length[x] < 3]];

(*generate side polygons - heights *)
SideComplex[pts_List, length_] := 
  Module[{topPts, botPts, sideRects, sidePts, sideNormals},
   topPts = pts;
   botPts = (2 length + 1 - #) & /@ topPts;
   sideRects = 
    Partition[
     RotateLeft[Flatten[{#, #} & /@ Range[Length[topPts]], 1]], 2];
   sidePts = {topPts[[#[[1]]]], botPts[[#[[1]]]], botPts[[#[[2]]]], 
       topPts[[#[[2]]]]} & /@ sideRects;
   Polygon@sidePts];

(* main code - it create top, bottom, and side polygons *)
To3DComplex[Polygon[list_], depth_: 10] := To3DComplex[list, depth]

To3DComplex[list_List, depth_: 10] /; (Depth[list] == 3) := 
 Module[{topPts, botPts, length, contours, sidePolys},
  topPts = {#[[1]], #[[2]], depth} & /@ list;
  botPts = Reverse[{#[[1]], #[[2]], 0} & /@ topPts];
  length = Length[list];
  contours = FindContourBreak[list];
  sidePolys = SideComplex[#, length] & /@ contours;
  GraphicsComplex[
   Join[topPts, botPts], {Polygon[Range[length]], 
    Polygon[Range[length + 1, 2 length] // Reverse], EdgeForm[], 
    sidePolys}]
  ]

To3DComplex[list_List, depth_: 10] := To3DComplex[#, depth] & /@ list

Here’s example:

divisions = 
  EntityValue[Entity["AdministrativeDivision", {_, "UnitedStates"}], 
   "Entities"];

project geoposition to mercator :

dat = (EntityValue[divisions, {"Population", "Polygon"}] /. 
      GeoPosition[x_] :> 
       GeoGridPosition[GeoPosition[x], "Mercator"]) /. 
    GeoGridPosition[x_, "Mercator"] :> x /. Quantity[x_, _] :> x;

rescale population for color function and depth:

pop = Rescale[(# - Min[#]) &@Log[dat[[All, 1]]] // N];

final result (I multiply 20 for depth):

poly = {ColorData["Rainbow"][#1], To3DComplex[#2, 20 #1]} & @@@ 
   Transpose[{pop, dat[[All, 2]]}];

Graphics3D[poly, PlotRange -> {{-60, -130}, {23, 60}, All}, 
 ImageSize -> 800, Boxed -> False]

ScreenShot2014-08-25at9.08.28AM.png
ScreenShot2014-08-25at9.08.28AM.png745×470 21.5 KB

9

That is very nice Jaebum!

On a related note, I’m trying getting up to speed with Entity[]. When going though this example, I was wondering how to delete members from an EntityList (e.g., Alaska and Hawaii). The following seems a bit cludgy to me

el = EntityList["Planet"]
DeleteCases[el, a_ /; MemberQ[a, "Venus"]]

I’ve glanced at the documentation, but haven’t found a way that would be more intuitive to beginners (i.e., my students). Is there a better way?

Perhaps this should be a separate thread?

10

you could do :

el = EntityList["Planet"];
DeleteCases[el, a_ /; MemberQ[a, "Venus", Infinity]]

divisions = 
 EntityValue[Entity["AdministrativeDivision", {_, "UnitedStates"}],   "Entities"];
DeleteCases[divisions,   a_ /; MemberQ[a, "Alaska" | "Hawaii", Infinity]]

for states you can also do:

divisions = 
 EntityValue[Entity["AdministrativeDivision", {Except["Alaska" | "Hawaii"],   "UnitedStates"}], "Entities"]
11

Yes. But, my goal here is different.

I am looking for something that is more intuitive for students and beginners. Yours are good examples for just-beyond-beginners who are learning about pattern restrictions. For beginners, wouldn’t DeleteEntities[entityList, {entityPattern, entityPattern}] be more initially more useful than writing it oneself?

12

DeleteCases do almost same thing you described:

DeleteCases[divisions, Entity[_, {"Alaska" | "Hawaii", _}]]

or

DeleteCases[el, Entity[_,"Venus"]]