如何生成包含最少n个点的不规则网格？

#1 楼

我看到MerseyViking建议使用四叉树。我将建议相同的内容，并且为了解释它，这是代码和示例。该代码是用R编写的，但是应该很容易地移植到Python中。 y方向，在每个级别上交替变化，直到不再需要分割为止。

由于目的是掩盖实际点的位置，因此在分割中引入一些随机性是很有用的。一种快速简单的方法是在分位数集上从50％处拆分一个小的随机量。以这种方式（a）分裂值极不可能与数据坐标重合，因此点将唯一地落入由分区创建的象限中，并且（b）点坐标将无法从四叉树精确地重建。 />
因为要在每个四叉树叶子中维持最小数量的节点数量k，所以我们实施了四叉树的受限形式。它将支持（1）将聚类点划分为每个在k和2 * k -1个元素之间的组，以及（2）映射象限。他们按班级。类标签加快了后处理（例如绘图）的速度，如下所示。该代码使用数字值作为ID。这在树中的深度达到52（使用双精度；如果使用无符号长整数，则最大深度为32）。对于较深的树（在任何应用中都不太可能，因为至少涉及R * 2 ^ 52点），所以id必须是字符串。

quadtree <- function(xy, k=1) {
  d = dim(xy)[2]
  quad <- function(xy, i, id=1) {
    if (length(xy) < 2*k*d) {
      rv = list(id=id, value=xy)
      class(rv) <- "quadtree.leaf"
    }
    else {
      q0 <- (1 + runif(1,min=-1/2,max=1/2)/dim(xy)[1])/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))
      class(rv) <- "quadtree"
    }
    return(rv)
  }
  quad(xy, 1)
}

请注意，此算法（以及大多数后处理算法）的递归分治设计意味着时间要求为O（m），RAM使用量为O（n），其中k是数字单元格数和m是点数。 n与m成比例，除以每个单元格的最小点n。这对于估计计算时间很有用。例如，如果将n = 10 ^ 6点划分为50-99点（k = 50）的单元格需要13秒，则m = 10 ^ 6/50 =20000。如果您想将其划分为5-9每个像元的点数（k = 5），m大10倍，因此计时时间增加到大约130秒。 （因为随着单元的变小，在它们的中心附近分解一组坐标的过程变得更快，所以实际的时间只有90秒。）要一直达到每个单元的k = 1点，大约需要六倍的时间。还是九分钟，我们可以预期代码实际上会比这快一点。

在继续之前，让我们生成一些有趣的不规则空间数据并创建其受限的四叉树（0.29秒的经过时间）：



下面是生成这些图的代码。它利用了k的多态性：例如，只要将R函数应用于points.quadtree对象，就会调用points。函数的强大功能非常简单，可以根据聚类标识符为点着色，这很明显：需要重复裁剪用于四叉树划分的阈值，但是相同的递归方法既简单又优雅。如果需要，可以使用变体构造象限的多边形表示形式。

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

作为另一个示例，我生成了1,000,000个点并将它们分为5-9个组。计时是91.7秒。

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!

作为与GIS交互的一个示例，让我们使用quadtree库将所有四叉树单元作为多边形shapefile写出。该代码模拟shapefiles的裁剪例程，但是这一次它必须生成单元的向量描述。这些作为数据帧输出，供lines.quadtree库使用。

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=".")

点本身可以​​使用shapefiles或通过导入（x，y ）坐标。

使用示例：

cell <- function(q, xylim, ...) {
  if (class(q)=="quadtree") f <- cell.quadtree else f <- cell.quadtree.leaf
  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]))
}

（在此处对read.shp使用任何所需的范围以放大到子区域或扩展映射到更大的区域；此代码默认为点的范围。）

仅此一项就足够了：这些多边形与原始点的空间连接将识别聚类。一旦确定，数据库“摘要”操作将生成每个单元格内点的摘要统计信息。

#2 楼

查看此算法是否为您的数据样本提供了足够的匿名性：


从规则网格开始
如果多边形的阈值小于阈值，则与相邻的交替合并（E，S，W ，N）顺时针旋转。
如果多边形小于阈值，则转到2，否则转到下一个多边形

，例如，最小阈值为3：

#3 楼

与Paulo有趣的解决方案类似，如何使用四叉树细分算法？

设置希望四叉树去的深度。您还可以为每个像元设置最小或最大点数，以便某些节点比其他节点更深/更小。冲洗并重复直到符合条件。

#4 楼

另一种方法是创建一个非常精细的网格，并使用max-p算法。 http://pysal.readthedocs.org/en/v1.7/library/region/maxp.html