The purpose of this clustering is to simplify the display of point symbols: when many are close together on the map, they will be replaced by a single symbol to indicate a group.
The requirements point to a need for a simple, adaptive solution: the point symbols can be updated and as the user zooms in, different symbols will appear at different places on the map (or screen) extent.
An excellent candidate clearly is a region quadtree.
There is a simpler method that will act like a region quadtree. It requires less coding, no advance creation of a data structure, but you pay a (small) price by needing to perform some on-the-fly calculations during zooming and panning. Just grid the map. Specifically, suppose there are n point symbols to be drawn within the map's current extent which has a length of dx and a height of dy. Relative to the map's origin the symbols need to be drawn at coordinates (x[i], y[i]), i = 1, 2, ..., n. Selecting a grid cellsize of c partitions the map into a grid of cells. The cell in which the location ( x, y ) belongs is in row j (y) = Floor[y / c] and column j (x) (counting from 0, with rows going from bottom to top and columns from left to right). You can consider any cell receiving two or more points to be a "cluster". The cluster symbol can be drawn at the cell's center, which has coordinates.(j + c/2, k + c/2).
This leads to the following solution, presented in pseudocode:
m = Floor(dy/c)+1
n = Floor(dx/c)+1
Dimension a[m,n] = 0
For each (x[i], y[i]) to be displayed:
Increment( a[ j(y[i]), j(x[i]) ] )
End for
For each (x[i], y[i]) to be displayed:
row = j(y[i])
col = j(x[i])
If a[row, col] > 1:
Draw a symbol for a cluster of k points at (c*(col+0.5), c*(row+0.5))
a[row, col] = 0
Else
Draw a point symbol at (x[i], y[i])
End if
End for
Clearly the algorithm's computing burden is O(# of points) in time and O(dx*dy/c^2) in storage. The trade-offs involved in selecting the cellsize c are:
c should be as large as possible: Storage is inversely proportional to c^2: small values of c mean large amounts of RAM. (Storage can be reduced to O(# of points) by using sparse arrays or dictionaries.)
c should be as large as possible: Two symbols (points or clusters) will never be any closer than c/2.
c should be as small as possible: every cluster symbol represents locations no further than c/sqrt(2) away from it.
c should be as small as possible: Large values of c tend to create many clusters and allow few individual points to appear.
Let's do a quick analysis of (4). As a point of departure, suppose the locations to be mapped occur uniformly at random and independently of one another. The number of cells is N (c) = (Floor(dx / c)+1) * (Floor(dy / c)+1), which--at least for larger values of c--is directly proportional to c^2. The distribution of cell counts will approximately follow a Poisson law with intensity lambda = n / N( c ) = n * c^2 / ( dx * dy ). The expected number of clusters equals
e (c) = n (1 - exp(-lambda)(1 + lambda)).
This becomes smaller as lambda shrinks to 0; that is, as the cellsize c gets smaller and smaller. The point of this analysis is that the preceding formula lets you anticipate how many clusters there might be, so you can select an initial value of c for which e (c) is below an acceptable value (while still being large enough to limit the RAM requirements). There is no closed-form solution, but a few Newton-Raphson steps will converge to one rapidly.
This approach is so dynamic--it's fast enough that the gridding and consequent clustering can be computed with each zoom and pan, and requires no precomputed data structures--that the desired "incremental modifications" as data are updated will happen automatically.