How can I generate irregular grid containing minimum n points?

Given a large (~ 1 million) sample of unevenly distributed points - is it possible to generate irregular grid (in size, but could also be irregular in shape if that is possible?) that will contain specified minimum amount of n points?

It is of less importance for me if genereted 'cells' of such grid contain exactly n number of points or at least n points.

I'm aware of solutions like genvecgrid in ArcGIS or Create Grid Layer in QGIS/mmgis however they will all create regular grids which will result in an output with empty cells (smaller problem - I could simply discard them) or cells with count of points less than n (bigger problem since I'd need a solution to aggregate those cells, probably using some tools from here?).

I've been googling around to no avail and am open to both commercial (ArcGIS & extensions) or free (Python, PostGIS, R) solutions.

• How "regular" does the grid need to be? I wonder if you can do some hierarchical clustering and then just cut the dendrogram to meet your needs (although this probably stretches what would be defined as a regular spatial configuration). The CrimeStat documentation has some good examples of this type of clustering. – Andy W Aug 10 '12 at 14:46
• Could you please explain exactly what you mean by an "irregular grid"? That sounds like an oxymoron :-). More to the point, what would be the purpose of this exercise? Note, too, that additional criteria or constraints are likely needed: after all, if you drew a square around all 1 million points, it could be considered as part of a grid and it would contain more than n of them. You probably wouldn't care for this trivial solution, though: but why not, exactly? – whuber Aug 10 '12 at 15:01
• @AndyW Thanks. Good idea & worth exploring. Will have a look. Size & shape of 'grid' is of secondary importance for me - priority (due to data privacy) is to 'hide' n features behind one – radek Aug 10 '12 at 16:18
• @whuber Thanks as well. I do agree - but wasn't sure how else I could name such partitioning. As mentioned above - my main motivation is data privacy. Having five point locations (which I cannot show on final map) I'd like to represent them by area covering them; and get mean/median/etc. value for that. I agree that it would be possible to draw one rectangle or convex hull representing them all - that would be the ultimate data privacy protection I guess? ;] However - it'd be more useful to represent it by shapes bounding, let's say 10 features. Then - I can still preserve spatial pattern. – radek Aug 10 '12 at 16:23
• IMO given your description I would utilize some type of interpolation and display a raster map (perhaps an adaptive bandwidth the size of your minimal N would be sufficient to smooth the data). As far as CrimeStat, the largest files I have used were around 100,000 cases I believe (and the clustering would certainly take along time). It is likely you could do some pre-generalization of your data to represent it as fewer cases and still get desirable results for whatever you want. It is a really simple program, I would suggest just taking the few minutes to try it out and see for yourself. – Andy W Aug 10 '12 at 16:47

I see MerseyViking has recommended a quadtree. I was going to suggest the same thing and in order to explain it, here's the code and an example. The code is written in R but ought to port easily to, say, Python.

The idea is remarkably simple: split the points approximately in half in the x-direction, then recursively split the two halves along the y-direction, alternating directions at each level, until no more splitting is desired.

Because the intent is to disguise actual point locations, it is useful to introduce some randomness into the splits. One fast simple way to do this is to split at a quantile set a small random amount away from 50%. In this fashion (a) the splitting values are highly unlikely to coincide with data coordinates, so that the points will fall uniquely into quadrants created by the partitioning, and (b) point coordinates will be impossible to reconstruct precisely from the quadtree.

Because the intention is to maintain a minimum quantity k of nodes within each quadtree leaf, we implement a restricted form of quadtree. It will support (1) clustering points into groups having between k and 2*k-1 elements each and (2) mapping the quadrants.

This R code creates a tree of nodes and terminal leaves, distinguishing them by class. The class labeling expedites post-processing such as plotting, shown below. The code uses numeric values for the ids. This works up to depths of 52 in the tree (using doubles; if unsigned long integers are used, the maximum depth is 32). For deeper trees (which are highly unlikely in any application, because at least k * 2^52 points would be involved), ids would have to be strings.

quadtree <- function(xy, k=1) {
d = dim(xy)
quad <- function(xy, i, id=1) {
if (length(xy) < 2*k*d) {
rv = list(id=id, value=xy)
}
else {
q0 <- (1 + runif(1,min=-1/2,max=1/2)/dim(xy))/2 # Random quantile near the median
x0 <- quantile(xy[,i], q0)
j <- i %% d + 1 # (Works for octrees, too...)
rv <- list(index=i, threshold=x0,
lower=quad(xy[xy[,i] <= x0, ], j, id*2),
upper=quad(xy[xy[,i] > x0, ], j, id*2+1))
}
return(rv)
}
}

Note that the recursive divide-and-conquer design of this algorithm (and, consequently, of most of the post-processing algorithms) means that the time requirement is O(m) and RAM usage is O(n) where m is the number of cells and n is the number of points. m is proportional to n divided by the minimum points per cell, k. This is useful for estimating computation times. For instance, if it takes 13 seconds to partition n=10^6 points into cells of 50-99 points (k=50), m = 10^6/50 = 20000. If you want instead to partition down to 5-9 points per cell (k=5), m is 10 times larger, so the timing goes up to about 130 seconds. (Because the process of splitting a set of coordinates around their middles gets faster as the cells get smaller, the actual timing was only 90 seconds.) To go all the way to k=1 point per cell, it will take about six times longer still, or nine minutes, and we can expect the code actually to be a little faster than that.

Before going further, let's generate some interesting irregularly spaced data and create their restricted quadtree (0.29 seconds elapsed time): Here's the code to produce these plots. It exploits R's polymorphism: points.quadtree will be called whenever the points function is applied to a quadtree object, for instance. The power of this is evident in the extreme simplicity of the function to color the points according to their cluster identifier:

points.quadtree <- function(q, ...) {
points(q\$lower, ...); points(q\$upper, ...)
}
points.quadtree.leaf <- function(q, ...) {
points(q\$value, col=hsv(q\$id), ...)
}

Plotting the grid itself is a little trickier because it requires repeated clipping of the thresholds used for the quadtree partitioning, but the same recursive approach is simple and elegant. Use a variant to construct polygonal representations of the quadrants if desired.

lines.quadtree <- function(q, xylim, ...) {
i <- q\$index
j <- 3 - q\$index
clip <- function(xylim.clip, i, upper) {
if (upper) xylim.clip[1, i] <- max(q\$threshold, xylim.clip[1,i]) else
xylim.clip[2,i] <- min(q\$threshold, xylim.clip[2,i])
xylim.clip
}
if(q\$threshold > xylim[1,i]) lines(q\$lower, clip(xylim, i, FALSE), ...)
if(q\$threshold < xylim[2,i]) lines(q\$upper, clip(xylim, i, TRUE), ...)
xlim <- xylim[, j]
xy <- cbind(c(q\$threshold, q\$threshold), xlim)
lines(xy[, order(i:j)],  ...)
}
lines.quadtree.leaf <- function(q, xylim, ...) {} # Nothing to do at leaves!

As another example, I generated 1,000,000 points and partitioned them into groups of 5-9 each. Timing was 91.7 seconds.

n <- 25000       # Points per cluster
n.centers <- 40  # Number of cluster centers
sd <- 1/2        # Standard deviation of each cluster
set.seed(17)
centers <- matrix(runif(n.centers*2, min=c(-90, 30), max=c(-75, 40)), ncol=2, byrow=TRUE)
xy <- matrix(apply(centers, 1, function(x) rnorm(n*2, mean=x, sd=sd)), ncol=2, byrow=TRUE)
k <- 5
system.time(qt <- quadtree(xy, k))
#
# Set up to map the full extent of the quadtree.
#
xylim <- cbind(x=c(min(xy[,1]), max(xy[,1])), y=c(min(xy[,2]), max(xy[,2])))
plot(xylim, type="n", xlab="x", ylab="y", main="Quadtree")
#
# This is all the code needed for the plot!
#
lines(qt, xylim, col="Gray")
points(qt, pch=".") As an example of how to interact with a GIS, let's write out all the quadtree cells as a polygon shapefile using the shapefiles library. The code emulates the clipping routines of lines.quadtree, but this time it has to generate vector descriptions of the cells. These are output as data frames for use with the shapefiles library.

cell <- function(q, xylim, ...) {
f(q, xylim, ...)
}
cell.quadtree <- function(q, xylim, ...) {
i <- q\$index
j <- 3 - q\$index
clip <- function(xylim.clip, i, upper) {
if (upper) xylim.clip[1, i] <- max(q\$threshold, xylim.clip[1,i]) else
xylim.clip[2,i] <- min(q\$threshold, xylim.clip[2,i])
xylim.clip
}
d <- data.frame(id=NULL, x=NULL, y=NULL)
if(q\$threshold > xylim[1,i]) d <- cell(q\$lower, clip(xylim, i, FALSE), ...)
if(q\$threshold < xylim[2,i]) d <- rbind(d, cell(q\$upper, clip(xylim, i, TRUE), ...))
d
}
cell.quadtree.leaf <- function(q, xylim) {
data.frame(id = q\$id,
x = c(xylim[1,1], xylim[2,1], xylim[2,1], xylim[1,1], xylim[1,1]),
y = c(xylim[1,2], xylim[1,2], xylim[2,2], xylim[2,2], xylim[1,2]))
}

The points themselves can be read directly using read.shp or by importing a data file of (x,y) coordinates.

Example of use:

qt <- quadtree(xy, k)
xylim <- cbind(x=c(min(xy[,1]), max(xy[,1])), y=c(min(xy[,2]), max(xy[,2])))
polys <- cell(qt, xylim)
polys.attr <- data.frame(id=unique(polys\$id))
library(shapefiles)
polys.shapefile <- convert.to.shapefile(polys, polys.attr, "id", 5)

(Use any desired extent for xylim here to window into a subregion or to expand the mapping to a larger region; this code defaults to the extent of the points.)

This alone is enough: a spatial join of these polygons to the original points will identify the clusters. Once identified, database "summarize" operations will generate summary statistics of the points within each cell.

• Wow! Fantastic. Will give it a shot with my data once back in the office =) – radek Aug 21 '12 at 15:48
• Top answer @whuber! +1 – MerseyViking Aug 21 '12 at 17:33
• (1) You can read shapefiles directly with (inter alia) the shapefiles package or else you can export (x,y) coordinates in ASCII text and read them with read.table. (2) I recommend writing out qt in two forms: first, as a point shapefile for xy where the id fields are included as cluster identifiers; second, where the line segments plotted by lines.quadtree are written out as a polyline shapefile (or where analogous processing writes the cells as a polygon shapefile). This is as simple as modifying lines.quadtree.leaf to output xylim as a rectangle. (See the edits.) – whuber Aug 21 '12 at 20:56
• @whubber Thanks a lot for an update. Everything worked smoothly. Well deserved +50, although now I think it deserves +500! – radek Aug 22 '12 at 10:16
• I suspect the calculated ids were not unique for some reason. Make these changes in the definition of quad: (1) initialize id=1; (2) change id/2 to id*2 in the lower= line; (3) make a similar change to id*2+1 in the upper= line. (I will edit my reply to reflect that.) That should also take care of the area calculation: depending on your GIS, all areas will be positive or all will be negative. If they're all negative, reverse the lists for x and y in cell.quadtree.leaf. – whuber Aug 22 '12 at 13:55

See if this algorithm gives enough anonymity for your data sample:

2. if polygon has less than threshold, merge with neighbor alternating (E, S, W, N) spiraling clockwise.
3. if polygon has less than threshold, go to 2, else go to next polygon

For example, if the minimum threshold is 3: • The devil is in the details: it seems this approach (or almost any agglomerative clustering algorithm) threatens to leave little "orphan" points scattered all over the place, which then cannot be processed. I'm not saying this approach is impossible, but I would maintain a healthy scepticism in the absence of an actual algorithm and example of its application to a realistic point dataset. – whuber Aug 21 '12 at 16:50
• Indeed this approach might be problematic. One application of this method that I was thinking about uses points as representations of residential buildings. I think this method would work well in more densely populated areas. However, there would be still cases when there is literally one or two buildings 'in the middle of nowhere' and it would take a lot of iterations & would result in really large areas to finally reach the minimum threshold. – radek Aug 21 '12 at 19:53

Similarly to Paulo's interesting solution, how about using a quad tree subdivision algorithm?

Set the depth you would like the quadtree to go to. You could also have a minimum or maximum number of points per cell so some nodes would be deeper/smaller than others.

Subdivide your world, discarding empty nodes. Rinse and repeat until the criteria are met.

• Thanks. What software would you recommend for that? – radek Aug 21 '12 at 10:41
• In principle this is a good idea. But how would empty nodes arise if you never allow less than a positive minimum number of points per cell? (There are many kinds of quadtrees, so the possibility of empty nodes indicates you have in mind one that is not adapted to the data, which raises concerns about its usefulness for the intended task.) – whuber Aug 21 '12 at 16:20
• I picture it like this: imagine a node has more than the maximum threshold of points in it, but they are clustered towards the top-left of the node. The node will be subdivided, but the bottom-right child node will be empty, so it can be pruned. – MerseyViking Aug 21 '12 at 17:31
• I see what you're doing (+1). The trick is to subdivide at a point determined by the coordinates (such as their median), thereby guaranteeing no empty cells. Otherwise, the quadtree is determined primarily by the space occupied by the points and not the points themselves; your approach then becomes an effective way of carrying out the generic idea proposed by @Paulo. – whuber Aug 21 '12 at 18:00

Another approach is to create a very fine grid, and use the max-p algorithm. http://pysal.readthedocs.org/en/v1.7/library/region/maxp.html