Sign up ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

I have a set of {x, y} coordinates, for example:

xmin = 0;
xmax = 100;
ymin = 42;
ymax = 76;
numPoints = 10^5;
pointList = Table[{RandomReal[{xmin, xmax}], RandomReal[{ymin, ymax}]}, {numPoints}];

I would like to break up this set of points into N equal-area squares in the most efficient manner possible, and without binning a point into more than one squares, such that we return a list like the following:

decomposedPointList = 
  {{...Points In Box 1...},{...Points In Box 2...},...,{...Points In Box N...};

Where the intersection of any two sublists of points is always zero.

Is there a nice way to do this?

share|improve this question

marked as duplicate by Kuba, Sjoerd C. de Vries, m_goldberg, Mr.Wizard Aug 14 '13 at 15:02

This question has been asked before and already has an answer. If those answers do not fully address your question, please ask a new question.

Did you really mean to write pointList=Table[{RandomReal[{xmin,xmax}],RandomReal[ymin,ymax]},{i,1,numPoints}‌​]; ? Should this not be pointList = Table[{RandomReal[{xmin, xmax}], RandomReal[{ymin, ymax}]}, {numPoints}]; so that the second entry is a point and not a list? also, your question is not clear, at least to me. – Nasser Aug 14 '13 at 3:36
@Nasser Thanks for that catch, I rewrote the line of code. – T.T. Aug 14 '13 at 3:38
@Nasser Basically I want to take my point like, and partition points into sublists based on these points falling into non-overlapping boxes (where we specify either the size of the boxes or the number of boxes). – T.T. Aug 14 '13 at 3:39
If I understand correctly this is a duplicate of Partitioning a 2D... – Kuba Aug 14 '13 at 4:50
@gpap I agree, I've posted my earlier than you have provided this link and I wasn't paying attention, why haven't you flagged this question as duplicate then? – Kuba Aug 14 '13 at 12:09

1 Answer 1

This may not scale well but perhaps can be a start.

f[x1_, x2_, y1_, y2_, size_, pt_] := Module[
  {nx, ny, tab, np, dist},
  nx = (x2 - x1)/size;
  ny = (y2 - y1)/size;
  tab = Flatten[
    Table[{x1, y1} + 0.5 {i size, j size}, {i, 1, 2 nx, 2}, {j, 1, 
      2 ny, 2}], 1];
  np = First@Nearest[tab, pt];
  dist = EuclideanDistance[np, pt];
  If[dist <= size/2, {np, size},
      Interval[{np[[1]] - 0.5  size, np[[1]] + 0.5 size}], pt[[1]]] &&

      Interval[{np[[2]] - 0.5  size, np[[2]] + 0.5 size}], 
      pt[[2]]], {np, size}, {"not in region of interest",""}]

Essentially partitions rectangular interest into sqaures of side 'size' and tabulates the centre. The final argument is the test point. There will be a unique nearest center. If the distance is size/2 then it is clearly in this square. The corner regions are dealt with intervals. Points outside region of interest will be coded.

I tested only on small test set 100 points of varying square sizes and seems to work but may not scale well. For example:

test = Thread[{RandomReal[{0, 100}, 100], RandomReal[{42, 76}, 100]}];
     "Square centred at `1` of size `2`:", #[[1, 1]], #[[1, 
       2]]], #[[2]]} & /@ 
  Tally[f[0, 100, 42, 76, 1, #] & /@ test] // Grid
share|improve this answer

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