¿Cómo puedo generar una cuadrícula irregular que contenga un mínimo de n puntos?

20

Dada una muestra grande (~ 1 millón) de puntos distribuidos de manera desigual, ¿es posible generar una cuadrícula irregular (en tamaño, pero también podría ser de forma irregular si es posible?) Que contendrá una cantidad mínima especificada de n puntos?

Para mí es menos importante si las 'células' generadas de dicha cuadrícula contienen exactamente n número de puntos o al menos n puntos.

Conozco soluciones como genvecgrid en ArcGIS o Create Grid Layer en QGIS / mmgis, sin embargo, todas crearán cuadrículas regulares que darán como resultado una salida con celdas vacías (problema más pequeño, simplemente podría descartarlas) o celdas con conteo de puntos menos de n (un problema mayor ya que necesitaría una solución para agregar esas celdas, ¿probablemente usando algunas herramientas de aquí ?).

He estado buscando en Google sin éxito y estoy abierto a soluciones comerciales (ArcGIS y extensiones) o gratuitas (Python, PostGIS, R).

radek
fuente
1
¿Qué tan "regular" debe ser la grilla? Me pregunto si puede hacer un agrupamiento jerárquico y luego simplemente cortar el dendrograma para satisfacer sus necesidades (aunque esto probablemente amplía lo que se definiría como una configuración espacial regular). La documentación de CrimeStat tiene algunos buenos ejemplos de este tipo de agrupación .
Andy W
55
¿Podría explicar exactamente qué quiere decir con una "cuadrícula irregular"? Eso suena como un oxímoron :-). Más concretamente, ¿cuál sería el propósito de este ejercicio? Tenga en cuenta también que es probable que se necesiten criterios o restricciones adicionales: después de todo, si dibujó un cuadrado alrededor de los 1 millón de puntos, podría considerarse como parte de una cuadrícula y contendría más de n de ellos. Sin embargo, probablemente no le interesaría esta solución trivial: pero ¿por qué no exactamente?
whuber
@AndyW Gracias. Buena idea y vale la pena explorar. Echaré un vistazo. Tamaño y forma de 'rejilla' es de importancia secundaria para mí - la prioridad (debido a la privacidad de los datos) es 'ocultar' n dispone detrás de uno
Radek
@whuber Gracias también. Estoy de acuerdo, pero no estaba seguro de qué otra forma podría nombrar esa partición. Como se mencionó anteriormente, mi principal motivación es la privacidad de los datos. Con cinco ubicaciones de puntos (que no puedo mostrar en el mapa final) me gustaría representarlas por área que las cubra; y obtener media / mediana / etc. valor para eso. Estoy de acuerdo en que sería posible dibujar un rectángulo o un casco convexo que los represente a todos: ¿sería la mejor protección de la privacidad de datos? ;] Sin embargo, sería más útil representarlo mediante delimitación de formas, digamos 10 características. Entonces, todavía puedo preservar el patrón espacial.
radek
1
OMI dada su descripción, utilizaría algún tipo de interpolación y mostraría un mapa ráster (quizás un ancho de banda adaptativo del tamaño de su N mínimo sería suficiente para suavizar los datos). En cuanto a CrimeStat, creo que los archivos más grandes que he usado fueron alrededor de 100,000 casos (y la agrupación ciertamente llevaría mucho tiempo). Es probable que pueda hacer una generalización previa de sus datos para representarlos como menos casos y aún así obtener resultados deseables para lo que desee. Es un programa realmente simple, sugeriría tomarse unos minutos para probarlo y verlo por usted mismo.
Andy W

Respuestas:

26

Veo que MerseyViking ha recomendado un quadtree . Iba a sugerir lo mismo y para explicarlo, aquí está el código y un ejemplo. El código está escrito Rpero debe portarse fácilmente a, por ejemplo, Python.

La idea es notablemente simple: divida los puntos aproximadamente a la mitad en la dirección x, luego divida recursivamente las dos mitades a lo largo de la dirección y, alternando las direcciones en cada nivel, hasta que no se desee más división.

Debido a que la intención es disfrazar ubicaciones de puntos reales, es útil introducir algo de aleatoriedad en las divisiones . Una manera rápida y sencilla de hacer esto es dividir en un conjunto de cuantiles una pequeña cantidad aleatoria del 50%. De esta manera (a) es poco probable que los valores de división coincidan con las coordenadas de datos, de modo que los puntos caerán únicamente en cuadrantes creados por la partición, y (b) las coordenadas de puntos serán imposibles de reconstruir con precisión a partir del quadtree.

Debido a que la intención es mantener una cantidad mínima kde nodos dentro de cada hoja de quadtree, implementamos una forma restringida de quadtree. Admitirá (1) puntos de agrupamiento en grupos que tienen entre k2 y k-1 elementos cada uno y (2) mapeo de los cuadrantes.

Este Rcódigo crea un árbol de nodos y hojas terminales, distinguiéndolos por clase. El etiquetado de la clase agiliza el procesamiento posterior, como el trazado, que se muestra a continuación. El código usa valores numéricos para los identificadores. Esto funciona hasta profundidades de 52 en el árbol (usando dobles; si se usan enteros largos sin signo, la profundidad máxima es 32). Para árboles más profundos (que son altamente improbables en cualquier aplicación, porque al menos k* 2 ^ 52 puntos estarían involucrados), los identificadores deberían ser cadenas.

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

Tenga en cuenta que el diseño recursivo de divide y vencerás de este algoritmo (y, en consecuencia, de la mayoría de los algoritmos de procesamiento posterior) significa que el requisito de tiempo es O (m) y el uso de RAM es O (n), donde mes el número de celdas y nes el número de puntos. mes proporcional a ndividido por los puntos mínimos por celda,k. Esto es útil para estimar los tiempos de cálculo. Por ejemplo, si se tarda 13 segundos en dividir n = 10 ^ 6 puntos en celdas de 50-99 puntos (k = 50), m = 10 ^ 6/50 = 20000. Si lo desea, en cambio, particione a 5-9 puntos por celda (k = 5), m es 10 veces más grande, por lo que el tiempo sube a unos 130 segundos. (Debido a que el proceso de dividir un conjunto de coordenadas alrededor de su centro se hace más rápido a medida que las celdas se hacen más pequeñas, el tiempo real fue de solo 90 segundos). Para llegar a k = 1 punto por celda, tomará aproximadamente seis veces más aún, o nueve minutos, y podemos esperar que el código sea un poco más rápido que eso.

Antes de continuar, generemos algunos datos interesantes espaciados irregularmente y creemos su quadtree restringido (tiempo transcurrido de 0.29 segundos):

Quadtree

Aquí está el código para producir estas parcelas. Explota Rel polimorfismo: points.quadtreese llamará siempre que la pointsfunción se aplique a un quadtreeobjeto, por ejemplo. El poder de esto es evidente en la extrema simplicidad de la función para colorear los puntos de acuerdo con su identificador de grupo:

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

Trazar la cuadrícula en sí es un poco más complicado porque requiere el recorte repetido de los umbrales utilizados para la partición de quadtree, pero el mismo enfoque recursivo es simple y elegante. Use una variante para construir representaciones poligonales de los cuadrantes si lo desea.

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 otro ejemplo, generé 1,000,000 de puntos y los dividí en grupos de 5-9 cada uno. El tiempo fue 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=".")

ingrese la descripción de la imagen aquí


Como un ejemplo de cómo interactuar con un SIG , escribamos todas las celdas del árbol cuádruple como un shapefilesarchivo de forma poligonal usando la biblioteca. El código emula las rutinas de recorte de lines.quadtree, pero esta vez tiene que generar descripciones vectoriales de las celdas. Estos se generan como marcos de datos para usar con la 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]))
}

Los puntos en sí pueden leerse directamente utilizando read.shpo importando un archivo de datos de coordenadas (x, y).

Ejemplo 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 cualquier extensión deseada xylimaquí para abrir una ventana en una subregión o para expandir el mapeo a una región más grande; este código se predetermina a la extensión de los puntos).

Esto por sí solo es suficiente: una unión espacial de estos polígonos a los puntos originales identificará los grupos. Una vez identificadas, las operaciones de "resumen" de la base de datos generarán estadísticas resumidas de los puntos dentro de cada celda.

whuber
fuente
¡Guauu! Fantástico. Lo
intentaré
44
¡La mejor respuesta @whuber! +1
MerseyViking
1
(1) Puede leer archivos de forma directamente con ( entre otros ) el shapefilespaquete o puede exportar coordenadas (x, y) en texto ASCII y leerlas con read.table. (2) Recomiendo escribir qten dos formas: primero, como un archivo de forma de punto xydonde los idcampos se incluyen como identificadores de clúster; segundo, donde los segmentos de línea trazados lines.quadtreese escriben como un archivo de forma de polilínea (o donde el procesamiento análogo escribe las celdas como un archivo de forma de polígono). Esto es tan simple como modificar lines.quadtree.leafa la salida xylimcomo un rectángulo. (Ver las ediciones.)
whuber
1
@whubber Muchas gracias por una actualización. Todo funcionó sin problemas. Bien merecido +50, aunque ahora creo que merece +500!
radek
1
Sospecho que los identificadores calculados no fueron únicos por alguna razón. Realice estos cambios en la definición de quad: (1) inicializar id=1; (2) cambiar id/2a id*2en la lower=línea; (3) hacer un cambio similar a id*2+1en la upper=línea. (Editaré mi respuesta para reflejar eso). Eso también debería ocuparse del cálculo del área: dependiendo de su SIG, todas las áreas serán positivas o todas serán negativas. Si todos son negativos, invierta las listas para xy ydentro cell.quadtree.leaf.
whuber
6

Vea si este algoritmo proporciona suficiente anonimato para su muestra de datos:

  1. comenzar con una cuadrícula regular
  2. Si el polígono tiene menos del umbral, combine con el vecino alternando (E, S, W, N) en espiral en el sentido de las agujas del reloj.
  3. si el polígono tiene menos del umbral, vaya al 2, de lo contrario vaya al siguiente polígono

Por ejemplo, si el umbral mínimo es 3:

algoritmo

Paulo Scardine
fuente
1
El diablo está en los detalles: parece que este enfoque (o casi cualquier algoritmo de agrupamiento aglomerativo) amenaza con dejar pequeños puntos "huérfanos" dispersos por todo el lugar, que luego no pueden procesarse. No digo que este enfoque sea imposible, pero mantendría un escepticismo saludable en ausencia de un algoritmo real y un ejemplo de su aplicación a un conjunto de datos de puntos realistas.
whuber
De hecho, este enfoque podría ser problemático. Una aplicación de este método en la que estaba pensando utiliza puntos como representaciones de edificios residenciales. Creo que este método funcionaría bien en áreas más densamente pobladas. Sin embargo, todavía habría casos en los que haya literalmente uno o dos edificios 'en el medio de la nada' y esto requeriría muchas iteraciones y daría como resultado áreas realmente grandes para finalmente alcanzar el umbral mínimo.
radek
5

De manera similar a la interesante solución de Paulo, ¿qué tal usar un algoritmo de subdivisión de árbol cuádruple?

Establezca la profundidad a la que desea que vaya el quadtree. También podría tener un número mínimo o máximo de puntos por celda para que algunos nodos sean más profundos / más pequeños que otros.

Subdivide su mundo, descartando nodos vacíos. Enjuague y repita hasta que se cumplan los criterios.

MerseyViking
fuente
Gracias. ¿Qué software recomendarías para eso?
radek
1
En principio esta es una buena idea. Pero, ¿cómo surgirían los nodos vacíos si nunca permite menos de un número mínimo positivo de puntos por celda? (Hay muchos tipos de quadtrees, por lo que la posibilidad de nodos vacío indica que tiene en mente una que no está adaptada a los datos, lo que plantea preocupaciones acerca de su utilidad para la tarea prevista.)
whuber
1
Me lo imagino así: imagine que un nodo tiene más que el umbral máximo de puntos, pero están agrupados hacia la parte superior izquierda del nodo. El nodo se subdividirá, pero el nodo secundario inferior derecho estará vacío, por lo que se puede podar.
MerseyViking
1
Veo lo que estás haciendo (+1). El truco es subdividir en un punto determinado por las coordenadas (como su mediana), garantizando así que no haya celdas vacías. De lo contrario, el quadtree está determinado principalmente por el espacio ocupado por los puntos y no por los puntos en sí; su enfoque se convierte en una forma efectiva de llevar a cabo la idea genérica propuesta por @Paulo.
whuber