Implementing hexagon binning in mathematica

Question

Hexagon bin plots are a useful way of visualising large datasets of bivariate data. Here are a few examples:

With bin frequency indicated by grey level...

..and by glyph size

There are packages for creating this kind of plot in both "R" and Python. Obviously, the idea is similar to DensityHistogram plots.

How would one go about generating hexagonal bins in Mathematica? Also, how would one control the size of a plotmarker based on the bin frequency?

Update

As a starting point I have tried to create a triangular grid of points:

vert1 = Table[{x, Sqrt[3] y}, {x, 0, 20}, {y, 0, 10}];
vert2 = Table[{1/2 x, Sqrt[3] /2 y}, {x, 1, 41, 2}, {y, 1, 21, 2}];

verttri = Flatten[Join[vert1, vert2], 1];

overlaying some data..

data = RandomReal[{0, 20}, {500, 2}];

ListPlot[{verttri, data}, AspectRatio -> 1]

next step might involve using Nearest:

nearbin = Nearest[verttri];

ListPlot[nearbin[#] & /@ data, AspectRatio -> 1]

This gives the location of vertices with nearby data points. Unfortunately, I can't see how to count those data points..

Answer 1

With the set-up you already have, you can do

nearbin = Nearest[Table[verttri[[i]] -> i, {i, Length@verttri}]];
counts = BinCounts[nearbin /@ data, {1, Length@verttri + 1, 1}];

which counts the number of data points nearest to each vertex. Then just draw the glyphs directly:

With[{maxCount = Max@counts}, 
 Graphics[
  Table[Disk[verttri[[i]], 0.5 Sqrt[counts[[i]]/maxCount]], {i, Length@verttri}],
  Axes -> True]]

The square root is so that the area of the glyphs, and the number of black pixels, corresponds to the number of data points in each bin. I used data = RandomVariate[MultinormalDistribution[{10, 10}, 7 IdentityMatrix[2]], 500] to get the following plot:

As Jens has commented already, though, this is a unnecessarily slow way of going about it. One ought to be able to directly compute the bin index from the coordinates of a data point without going through Nearest. This way was easy to implement and works fine for a 500-point dataset though.

---

Update: Here's an approach that doesn't require you to set up a background grid in advance. We'll directly find the nearest grid vertex for each data point and then tally them up. To do so, we'll break the hexagonal grid into rectangular tiles of size $1\times\sqrt3$. As it turns out, when you're in say the $[0,1]\times[0,\sqrt3]$ tile, your nearest grid vertex can only be one of the five vertices in the tile, $(0,0)$, $(1,0)$, $(1/2,\sqrt3/2)$, $(0,\sqrt3)$, and $(1,\sqrt3)$. We could work out the conditions explicitly, but let's just let Nearest do the work:

tileContaining[{x_, y_}] := {Floor[x], Sqrt[3] Floor[y/Sqrt[3]]};
nearestWithinTile = Nearest[{{0, 0}, {1, 0}, {1/2, Sqrt[3]/2}, {0, Sqrt[3]}, {1, Sqrt[3]}}];
nearest[point_] := Module[{tile, relative},
   tile = tileContaining[point];
   relative = point - tile;
   tile + First@nearestWithinTile[relative]];

The point is that a NearestFunction over just five points ought to be extremely cheap to evaluate—certainly much cheaper than your NearestFunction over the several hundred points in verttri. Then we just have to apply nearest on all the data points and tally the results.

tally = Tally[nearest /@ data];
With[{maxTally = Max[Last /@ tally]}, 
 Graphics[
  Disk[#[[1]], 1/2 Sqrt[#[[2]]/maxTally]] & /@ tally, 
  Axes -> True, AxesOrigin -> {0, 0}]]

Answer 2

As this answer got much more attention than I thought at first, I felt compelled to pack the answer as a function.

The following function draws three different hex-histogram representations of your data.

Code stolen liberally from the question and from Rahul's answer.

f[data_, cs_, ptype_] :=
 Module[{hc, vh, nearbin, counts, trr, maxCount},
  hc = Flatten[Table[{i, j}, (*hexagons positions*)
     {p, 0, 1},
     {j, 3 p/2 cs + Min[data[[All, 2]]], Max[data[[All, 2]]], 3 cs},
     {i, p Sqrt[3]/2 cs + Min[data[[All, 1]]], Max[data[[All, 1]]], 
      Sqrt[3] cs}], 2];

  vh = cs Vertices[Hexagon];
  nearbin = Nearest[Table[hc[[i]] -> i, {i, Length@hc}]];
  counts = BinCounts[nearbin /@ data, {1, Length@hc + 1, 1}];
  trr[v_, tr_] := Translate[Rotate[Polygon[v], Pi/2], tr];
  maxCount = Max@counts;

  Graphics[Table[
    Switch[ptype,
     1, {Opacity[counts[[n]]/maxCount], trr[vh, hc[[n]]]},
     2, trr[counts[[n]]/maxCount vh, hc[[n]]],
     3, trr[Sqrt[counts[[n]]/maxCount] vh, hc[[n]]]],
    {n, Length@hc}],
   Axes -> True]
  ]

<< Polytopes`
cs = 1/2;(*cell size*)
data = RandomVariate[MultinormalDistribution[{10, 10}, 7 IdentityMatrix[2]], 500];
GraphicsRow@Table[f[data, cs, i], {i, 3}]

The old answer can be found in the post's edit history - removed

Answer 3

Here is a standalone function incorporating belisarius' standalone framework and Rahul Narain's fast use of NearestFunction. Full credit to these guys!

 Needs["Polytopes`"];
 honeybin[data_, cs_, ptype_, {{xmin_, xmax_}, {ymin_, ymax_}}] :=

 (*"cs_" is the width of a bin,i.e.the distance between bin centers*)

  Module[{tileContaining, nearestWithinTile, nearest, tally, vh, hexr,
    trr},

   tileContaining[{x_, y_}] := {Floor[x], Sqrt[3] Floor[y/Sqrt[3]]};
   nearestWithinTile = 
     Nearest[{{0, 0}, {1, 0}, {1/2, Sqrt[3]/2}, {0, Sqrt[3]}, {1, Sqrt[3]}}];
   nearest[point_] := Module[{tile, relative}, tile = tileContaining[point];
   relative = point - tile;
     tile + First@nearestWithinTile[relative]];

   vh = cs 1.05 Vertices[Hexagon]/Sqrt[3];
   trr[v_, tr_] := Translate[Rotate[Polygon[v], Pi/2], tr];
   tally = Tally[cs (nearest /@ (data/cs))];

  With[{maxTally = Max[Last /@ tally]},  
   Graphics[Table[
     Switch[ptype,
      1, {Blend[{Lighter[Red, 0.99], Darker[Red, 0.6]},
         Sqrt[Last@tally[[n]]/maxTally]], trr[vh, First@tally[[n]]]},
      2, trr[Last@tally[[n]]/maxTally vh, First@tally[[n]]], 
      3, trr[Sqrt[Last@tally[[n]]/maxTally] vh, First@tally[[n]]]], 
         {n, Length@tally}], 
    Frame -> True, 
    PlotRange -> {{xmin, xmax}, {ymin, ymax}}, 
    PlotRangeClipping -> True]]]

note there is a fudge factor (1.05) to get rid of whitespace between the cells.

With a little Colour...

data = RandomVariate[MultinormalDistribution[{10, 10}, 7 IdentityMatrix[2]], 1000];
honeybin[data, 1, 1, {{1, 15}, {3, 20}}]