Efficiently counting the number of data points that fall within each square of a grid

Question

The following code works, but is inefficient. I would like to find a better implementation.

The problem is the following. I have a list of points, representing a matrix of 2 columns and N rows. Each point represents an event occurring on a map. I create a tessellation for the given map, creating a grid of x-by-y squares I want to calculate the number of events that fall into each square.

The result should be an x-by-y matrix where an element gives the number of events falling inside a given square. I will use this to do some more point pattern analysis, etc.

grid is created as a list of polygons, that correspond to the squares of the tessellation. It is created as follows, with x0 and y0 defining the “origin corner” of the map from which to start the tessellation:

grid = 
  Table[
    Polygon[
      {{x0 + (i - 1)*dsq, y0 + (j - 1)*dsq}, 
       {x0 + i*dsq, y0 + (j - 1)*dsq}, 
       {x0 + i*dsq, y0 + j*dsq}, 
       {x0 + (i - 1)*dsq, y0 + j*dsq}, 
       {x0 + (i - 1)*dsq, y0 + (j - 1)*dsq}}], 
    {i, 1, xGrid}, {j, 1, yGrid}];

{{Polygon[{{41.645, -87.884}, {41.6828, -87.884}, {41.6828, -87.8463}, 
           {41.645, -87.8463}, {41.645, -87.884}}], 
  Polygon[{{41.645, -87.8463}, {41.6828, -87.8463}, {41.6828, -87.8085}, 
           {41.645, -87.8085}, {41.645, -87.8463}}], 
  Polygon[{{41.645, -87.8085}, {41.6828, -87.8085}, {41.6828, -87.7708}, 
           {41.645, -87.7708}, {41.645, -87.8085}}], 
  ...}}

I initialize the matrix to store the results (resMat)

xGrid = 10;
yGrid = 10;
resMat = Table[0, {x, 1, xGrid}, {y, 1, yGrid}] ;

This is the expression I use to go through each square and calculate how many points are in it:

Table[
  If[RegionMember[grid[[x, y]], points[[n, All]]], 
    resMat[[x, y]] = resMat [[x, y]] + 1], 
  {n, 1, Length[points]}, {x, 1, xGrid}, {y, 1, yGrid}] ;

I use the following to quickly visualize the results:

ArrayPlot[Reverse[resMat], ColorFunction -> ColorData["SolarColors"]]

The computation of resMat works, but is really slow and inefficient. The RegionMember[…] is the bottleneck: for a grid size of 100-by-100, the analysis of a single point can take 10-20 sec. My point list contains hundreds of points, making the computation with the current code too time consuming.

I could not come up with a better implementation so far. Therefore, I will be glad to get some help.

Answer 1

Not as short and quick as Michael E2's solution, but it collects the indices of the events for each grid cell seperately. That might get useful for later analysis.

n = 10;    
xlist = Subdivide[0., 1., n];
ylist = Subdivide[0., 1., n];
centers = 
  Partition[
   Tuples[{MovingAverage[xlist, 2], MovingAverage[ylist, 2]}], n];

pts = RandomReal[{0, 1}, {1000000, 2}];
data = Partition[Nearest[
      pts -> Automatic, Flatten[centers, 1], {\[Infinity], 0.5/n}, 
      DistanceFunction -> ChessboardDistance],
     n
     ]; // AbsoluteTiming // First

0.210811

Now data[[i,j]] contains the indices of all pts that lie in the cell with center centers[[i,j]]. Here a visualization:

i = 4;
j = 7;
Graphics[{
  Darker@Green,
  EdgeForm[{Darker@Darker@Green, Thickness[0.025]}],
  Rectangle[{xlist[[i]], ylist[[j]]}, {xlist[[i]] + 1/n, ylist[[j]] + 1/n}],
  Red,
  PointSize[0.0001],
  Point[pts[[data[[i, j]]]]]
  },
 PlotRange -> {{0, 1}, {0, 1}},
 Frame -> True,
 GridLines -> {xlist, ylist}
 ]

You can obtain the diagram you actually asked for by

MatrixPlot[Map[Length, data, {2}]]

Remark

Alternatively to using Nearest, you could also use

data1 = BinLists[pts, {0, 1, 1/n}, {0, 1, 1/n}];

which is almost exactly as fast (and probably implemented precisely like the method above). Its disadvantage is that you don't get the indices but the points themselves, i.e., we have the relation

Sort[data1[[i, j]]] == Sort[pts[[data[[i, j]]]]]

True

Not getting the indices might make it needlessly complicated to gather metadata associated to the events contained in a cell.

Answer 2

Perhaps BinCounts?

Example from the documentation with 1000 points:

Count random pairs in bins of width 0.25 in both dimensions:

BinCounts[RandomReal[{-1, 1}, {1000, 2}], {-1, 1, .25}, {-1, 1, .25}]
(*
{{8, 14, 21, 14, 15, 14, 18, 17}, {16, 15, 17, 19, 13, 12, 13, 12}, {19, 13, 
  33, 16, 9, 17, 14, 13}, {14, 18, 14, 8, 18, 11, 11, 9}, {24, 20, 17, 13, 21,
   16, 16, 20}, {19, 23, 17, 14, 13, 14, 12, 21}, {17, 15, 14, 11, 14, 14, 18,
   14}, {20, 15, 17, 10, 15, 19, 16, 16}}
*)

10^6 points in a 100 x 100 grid:

BinCounts[RandomReal[{-1, 1}, {1000000, 2}], {-1, 1, .02}, {-1, 1, .02}] // 
  Dimensions // AbsoluteTiming
(*{0.09895, {100, 100}}  *)

Answer 3

If you only want counts on a regular grid, you can do a little better than BinCounts. Using @Michael's example:

SeedRandom[1]
data = RandomReal[{-1, 1}, {10^6, 2}];

r1 = BinCounts[data, {-1, 1, .02}, {-1, 1, .02}]; // RepeatedTiming

r2 = Partition[
    Values @ KeySort @ Counts[Floor[50 data] . {100, 1}],
    100
]; //RepeatedTiming

r1 === r2

{0.048, Null}

{0.032, Null}

True

Answer 4

I liked the idea Henrik Schumacher, with the Nearest approach, because I need to collect the indeces for the x and y axes, rather than simply count the number of points in each bin. Here below is my implementation:

x0 = 0;
y0 = 0;
xyPoints = {RandomReal[{0, 100}, 100], RandomReal[{0, 100}, 100]} // 
   Transpose;
tessSize = 10;
xGrid = 100/tessSize;
yGrid = 100/tessSize;

tessSize is the size of tesselation desired and {x0, y0} is the "origin" of the frame/image. Rather then defining a list of polygons, I define the steps of the tesselation in the x an y axis, as follow:

xSteps = Table[x0 + (i - 1)*tessSize, {i, 1, xGrid}]
ySteps = Table[y0 + (j - 1)*tessSize, {j, 1, yGrid}]

{0, 10, 20, 30, 40, 50, 60, 70, 80, 90}

{0, 10, 20, 30, 40, 50, 60, 70, 80, 90}

The set of points (xyPoints) is given in the form of a two-column matrix: column one for x and two for y coordinate. For a point we can get the position to the nearest tesselation point, xSteps, for example:

Position[xSteps,_?(#==Nearest[xSteps,xyPoints[[ii,1]]][[1]]&)]

Any given point will follow between two points. We can then discriminate to which tesselation square to asign by calculating the distance to the two closest tesslation points (Nearest). If distance > 0 (Boole[ ... ] == True), you found the right x coordinate of the grid - square. If distance <= 0 (Boole[ ... ] == False) then the right x coordinate is the previous position. To make sure that the list of indeces start at 1 and finish at Length[xSteps], we add: (Boole[ ... ]-1) + Postion[ ... ] Here is the core of the algorithm, the final implementation to find the indeces of tesselation for x an y, stored in a matrix

resMat = Transpose[{
   Flatten[
    {Table[Flatten[
       (Boole[ Positive[ xyPoints[[ii, 1]] - Nearest[xSteps, xyPoints[[ii, 1]]]]] - 1) + 
        Position[ xSteps, _?(# == Nearest[xSteps, xyPoints[[ii, 1]]][[1]] &)]], 
        {ii, 1, Length[xyPoints]}] }],

   Flatten[
    {Table[Flatten[
       (Boole[ Positive[ xyPoints[[ii, 2]] - Nearest[ySteps, xyPoints[[ii, 2]]]]] - 1) + 
        Position[ ySteps, _?(# ==  Nearest[ySteps, xyPoints[[ii, 2]]][[1]] &)]], 
        {ii, 1, Length[xyPoints]}] }]
   }]

{{1, 6}, {10, 1}, {5, 2}, {1, 10}, {10, 7}, {2, 7}, {10, 4}, {1, 8}, ... <<<90>>> ... {3, 4}, {5, 6}}

Lastly, a simple graph to see what is doing:

Show[
 ListPlot[ xyPoints, PlotStyle -> PointSize[Medium], AspectRatio -> 1, 
  Axes -> False],
 ListLinePlot[ Transpose[{Transpose[{xSteps, Table[0, Length[xSteps]]}], {xSteps,Table[100, Length[xSteps]]}\\Transpose}], 
  PlotStyle -> { {Red, Dashing[{0.005, 0.01}]}}], 
 ListLinePlot[ Transpose[{Transpose[{Table[0, Length[ySteps]], ySteps}], {Table[100, Length[ySteps]], ySteps}\\Transpose}], 
  PlotStyle -> { {Red, Dashing[{0.005, 0.01}]}}]
 ]

Thank you all for the help ;)