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

20

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

radek
fonte
1
Quão "regular" a grade precisa ser? Gostaria de saber se você pode fazer alguns agrupamentos hierárquicos e depois cortar o dendograma para atender às suas necessidades (embora isso provavelmente estenda o que seria definido como uma configuração espacial regular). A documentação do CrimeStat tem alguns bons exemplos desse tipo de cluster .
217 Andy W no dia
5
Você poderia, por favor, explicar exatamente o que você quer dizer com "grade irregular"? Isso soa como um oxímoro :-). Mais precisamente, qual seria o objetivo deste exercício? Observe também que provavelmente são necessários critérios ou restrições adicionais: afinal, se você desenhar um quadrado em torno de 1 milhão de pontos, ele poderá ser considerado como parte de uma grade e conterá mais de n deles. Você provavelmente não se importaria com esta solução trivial: mas por que não, exatamente?
whuber
@AndyW Thanks. Boa ideia e vale a pena explorar. Vai dar uma olhada. Tamanho e forma de 'grid' é de importância secundária para mim - prioridade (devido a privacidade de dados) é a 'esconder' n apresenta atrás de um
radek
@whuber Obrigado também. Eu concordo - mas não tinha certeza de como mais poderia nomear esse particionamento. Como mencionado acima - minha principal motivação é a privacidade dos dados. Tendo cinco localizações de pontos (que não posso mostrar no mapa final), gostaria de representá-las por área que as cobre; e obtenha média / mediana / etc. valor para isso. Concordo que seria possível desenhar um retângulo ou casco convexo representando todos eles - essa seria a melhor proteção de privacidade de dados, eu acho? ;] No entanto - seria mais útil representá-lo por formas delimitadas, digamos 10 recursos. Então - ainda posso preservar o padrão espacial.
Radek
1
Na IMO, dada sua descrição, eu usaria algum tipo de interpolação e exibia um mapa raster (talvez uma largura de banda adaptável do tamanho do seu N mínimo fosse suficiente para suavizar os dados). Quanto ao CrimeStat, os maiores arquivos que usei foram cerca de 100.000 casos, acredito (e o cluster certamente levaria tempo). É provável que você possa fazer uma pré-generalização de seus dados para representá-los como menos casos e ainda obter resultados desejáveis ​​para o que quiser. É um programa realmente simples, sugiro apenas alguns minutos para testá-lo e ver por si mesmo.
Andy W

Respostas:

26

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

Quadtree

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

insira a descrição da imagem aqui


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.

whuber
fonte
Uau! Fantástico. Vai dar-lhe um tiro com meus dados, uma vez de volta ao escritório =)
radek
4
Melhor resposta @whuber! 1)
MerseyViking
1
(1) Você pode ler os shapefiles diretamente com ( entre outros ) o shapefilespacote ou pode exportar (x, y) coordenadas no texto ASCII e lê-las com read.table. (2) eu recomendo escrever de qtduas formas: primeiro, como um shapefile de ponto para xyonde os idcampos são incluídos como identificadores de cluster; segundo, onde os segmentos de linha plotados por lines.quadtreesão gravados como um arquivo de forma polilinha (ou onde o processamento análogo grava as células como um arquivo de forma poligonal). Isso é tão simples quanto modificar lines.quadtree.leafa saída xylimcomo um retângulo. (Veja as edições.)
whuber
1
@whubber Muito obrigado por uma atualização. Tudo funcionou sem problemas. Bem merecido +50, embora agora eu acho que merece +500!
Radek
1
Suspeito que os IDs calculados não tenham sido exclusivos por algum motivo. Faça essas alterações na definição de quad: (1) inicializar id=1; (2) mude id/2para id*2na lower=linha; (3) faça uma alteração semelhante à id*2+1da upper=linha. (Vou editar minha resposta para refletir isso.) Isso também deve cuidar do cálculo da área: dependendo do seu GIS, todas as áreas serão positivas ou todas serão negativas. Se todas forem negativas, inverta as listas para xe para ydentro cell.quadtree.leaf.
whuber
6

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:

algoritmo

Paulo Scardine
fonte
1
O diabo está nos detalhes: parece que essa abordagem (ou quase qualquer algoritmo de agrupamento aglomerativo) ameaça deixar pequenos pontos "órfãos" espalhados por todo o lugar, que não podem ser processados. Não estou dizendo que essa abordagem seja impossível, mas manteria um ceticismo saudável na ausência de um algoritmo real e exemplo de sua aplicação a um conjunto de dados de pontos realistas.
whuber
De fato, essa abordagem pode ser problemática. Uma aplicação desse método em que eu estava pensando usa pontos como representações de edifícios residenciais. Eu acho que esse método funcionaria bem em áreas mais densamente povoadas. No entanto, ainda haveria casos em que houvesse literalmente um ou dois edifícios 'no meio do nada' e seriam necessárias muitas iterações e resultariam em áreas realmente grandes para finalmente atingir o limite mínimo.
Radek
5

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.

MerseyViking
fonte
Obrigado. Qual software você recomendaria para isso?
Radek
1
Em princípio, é uma boa ideia. Mas como surgiriam nós vazios se você nunca permitir menos que um número mínimo positivo de pontos por célula? (Existem muitos tipos de quadtrees, então a possibilidade de nós vazios indica que você tem em mente que não está adaptado aos dados, o que levanta preocupações sobre a sua utilidade para a tarefa pretendida.)
whuber
1
Eu imagino o seguinte: imagine que um nó tenha mais do que o limite máximo de pontos, mas eles estão agrupados na parte superior esquerda do nó. O nó será subdividido, mas o nó filho no canto inferior direito estará vazio, para que possa ser removido.
usar o seguinte código
1
Eu vejo o que você está fazendo (+1). O truque é subdividir em um ponto determinado pelas coordenadas (como a mediana), garantindo assim nenhuma célula vazia. Caso contrário, o quadtree é determinado principalmente pelo espaço ocupado pelos pontos e não pelos pontos; sua abordagem se torna uma maneira eficaz de levar a cabo a idéia genérica proposta por @Paulo.
whuber