Como posso gerar uma grade irregular contendo no mínimo n pontos?

Dada uma amostra grande (~ 1 milhão) de pontos desigualmente distribuídos - é possível gerar uma grade irregular (em tamanho, mas também pode ter uma forma irregular se for possível?) Que conterá uma quantidade mínima especificada de n pontos?

É menos importante para mim se as 'células' geradas dessa grade contenham exatamente n número de pontos ou pelo menos n pontos.

Estou ciente de soluções como genvecgrid no ArcGIS ou Create Grid Layer no QGIS / mmgis, no entanto, todas elas criarão grades regulares que resultarão em uma saída com células vazias (problema menor - eu poderia simplesmente descartá-las) ou células com contagem de pontos menor que n (problema maior, pois eu precisaria de uma solução para agregar essas células, provavelmente usando algumas ferramentas daqui ?).

Estou pesquisando sem sucesso e estou aberto a soluções comerciais (ArcGIS e extensões) ou gratuitas (Python, PostGIS, R).

Vejo que MerseyViking recomendou um quadtree . Eu sugeriria a mesma coisa e, para explicar, aqui está o código e um exemplo. O código está escrito, Rmas deve portar facilmente para, por exemplo, Python.

A idéia é notavelmente simples: divida os pontos aproximadamente ao meio na direção x, depois divida recursivamente as duas metades ao longo da direção y, alternando as direções em cada nível, até que não seja mais necessária a divisão.

Como a intenção é disfarçar a localização real dos pontos, é útil introduzir alguma aleatoriedade nas divisões . Uma maneira rápida e simples de fazer isso é dividir em um conjunto de quantis uma pequena quantidade aleatória de 50%. Dessa maneira (a) é altamente improvável que os valores de divisão coincidam com as coordenadas dos dados, de modo que os pontos caem exclusivamente nos quadrantes criados pela partição e (b) será impossível reconstruir precisamente as coordenadas dos pontos a partir da quadtree.

Como a intenção é manter uma quantidade mínima kde nós dentro de cada folha de quadtree, implementamos uma forma restrita de quadtree. Ele suportará (1) pontos de agrupamento em grupos que possuem entre ke 2 * k-1 elementos cada e (2) mapear os quadrantes.

Esse Rcódigo cria uma árvore de nós e folhas de terminal, distinguindo-os por classe. A rotulagem da classe agiliza o pós-processamento, como plotagem, mostrada abaixo. O código usa valores numéricos para os IDs. Isso funciona até profundidades de 52 na árvore (usando dobras; se números inteiros longos não assinados forem usados, a profundidade máxima será 32). Para árvores mais profundas (que são altamente improváveis ​​em qualquer aplicação, porque pelo menos k* 2 ^ 52 pontos estariam envolvidos), os IDs teriam que ser cadeias de caracteres.

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)
}

Observe que o design recursivo de dividir e conquistar deste algoritmo (e, conseqüentemente, da maioria dos algoritmos de pós-processamento) significa que o tempo necessário é O (m) e o uso de RAM é O (n) onde mé o número de células e né o número de pontos. mé proporcional ao ndividido pelos pontos mínimos por célula,k. Isso é útil para estimar os tempos de computação. Por exemplo, se levar 13 segundos para particionar n = 10 ^ 6 pontos em células de 50-99 pontos (k = 50), m = 10 ^ 6/50 = 20000. Se você desejar particionar para 5-9 pontos por célula (k = 5), m é 10 vezes maior, então o tempo sobe para cerca de 130 segundos. (Como o processo de dividir um conjunto de coordenadas em torno de seus médios fica mais rápido à medida que as células diminuem, o tempo real era de apenas 90 segundos.) Para percorrer todo o caminho de k = 1 ponto por célula, levará cerca de seis vezes mais ainda, ou nove minutos, e podemos esperar que o código seja realmente um pouco mais rápido que isso.

Antes de prosseguir, vamos gerar dados interessantes com espaçamento irregular e criar sua quadtree restrita (tempo decorrido de 0,29 segundos):

Aqui está o código para produzir esses gráficos. Ele explora Ro polimorfismo: points.quadtreeserá chamado sempre que a pointsfunção for aplicada a um quadtreeobjeto, por exemplo. O poder disso é evidente na extrema simplicidade da função de colorir os pontos de acordo com o identificador de cluster:

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

A plotagem da grade em si é um pouco mais complicada, pois requer recortes repetidos dos limites usados ​​para o particionamento de quadtree, mas a mesma abordagem recursiva é simples e elegante. Use uma variante para construir representações poligonais dos quadrantes, se desejado.

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!

Como outro exemplo, gerei 1.000.000 de pontos e os particionei em grupos de 5 a 9 cada. O tempo foi de 91,7 segundos.

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

Como um exemplo de como interagir com um GIS , vamos escrever todas as células quadtree como um arquivo de forma de polígono usando a shapefilesbiblioteca. O código emula as rotinas de recorte de lines.quadtree, mas desta vez ele precisa gerar descrições vetoriais das células. Eles são produzidos como quadros de dados para uso com a shapefilesbiblioteca.

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]))
}

Os pontos em si podem ser lidos diretamente usando read.shpou importando um arquivo de dados de coordenadas (x, y).

Exemplo de uso:

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)
write.shapefile(polys.shapefile, "f:/temp/quadtree", arcgis=TRUE)

(Use qualquer extensão desejada para exibir xylimaqui em uma sub-região ou expandir o mapeamento para uma região maior; esse código é padronizado na extensão dos pontos.)

Isso por si só é suficiente: uma junção espacial desses polígonos aos pontos originais identificará os agrupamentos. Uma vez identificadas, as operações de "resumo" do banco de dados geram estatísticas resumidas dos pontos dentro de cada célula.

Veja se este algoritmo fornece anonimato suficiente para sua amostra de dados:

1. comece com uma grade regular
2. se o polígono tiver menos que o limite, faça a mesclagem com o vizinho alternando (E, S, W, N) em espiral no sentido horário.
3. se o polígono tiver menos que o limite, vá para 2, caso contrário, vá para o próximo polígono

Por exemplo, se o limite mínimo for 3:

Da mesma forma que a interessante solução de Paulo, que tal usar um algoritmo de subdivisão de quatro árvores?

Defina a profundidade que você deseja que o quadtree vá. Você também pode ter um número mínimo ou máximo de pontos por célula, para que alguns nós sejam mais profundos / menores que outros.

Subdivida seu mundo, descartando nós vazios. Enxágue e repita até que os critérios sejam atendidos.

Outra abordagem é criar uma grade muito fina e usar o algoritmo max-p. http://pysal.readthedocs.org/en/v1.7/library/region/maxp.html