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

18

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 tener una forma irregular si eso es posible?) que contendrá una cantidad mínima especificada de n puntos?

Es de menor importancia para mí si las 'celdas' generadas de dicha cuadrícula contienen exactamente n número de puntos o al menos n puntos.

Tengo conocimiento de soluciones como genvecgrid en ArcGIS o Create Grid Layer en QGIS / mmgis, sin embargo, todos crearán grillas regulares que resultarán en una salida con celdas vacías (problema menor: simplemente podría descartarlos) o celdas con un conteo de puntos menor que n (problema mayor, ya que necesito una solución para agregar esas celdas, probablemente usando algunas herramientas de here ?).

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

    
pregunta radek 10.08.2012 - 14:56

4 respuestas

25

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 en R , pero se debe portar fácilmente a, digamos, Python.

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

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

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

Este código R crea un árbol de nodos y hojas de terminal, distinguiéndolos por clase. El etiquetado de clase acelera el procesamiento posterior, como el trazado, que se muestra a continuación. El código utiliza valores numéricos para los identificadores. Esto funciona hasta una profundidad 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, ya que estaría involucrado al menos k * 2 ^ 52 puntos), 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 división y conquista de este algoritmo (y, en consecuencia, de la mayoría de los algoritmos de posprocesamiento) significa que el requisito de tiempo es O (m) y el uso de RAM es O (n) donde m es el número de celdas y n es el número de puntos. m es proporcional a n dividido 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 particionar n = 10 ^ 6 puntos en celdas de 50-99 puntos (k = 50), m = 10 ^ 6/50 = 20000. Si, en cambio, desea particionar hasta 5-9 puntos por celda (k = 5), m es 10 veces más grande, por lo que el tiempo aumenta hasta aproximadamente 130 segundos. (Debido a que el proceso de dividir un conjunto de coordenadas alrededor de sus medianas partes se acelera 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. todavía, o nueve minutos, y podemos esperar que el código sea un poco más rápido que eso.

Antes de continuar, vamos a generar algunos datos interesantes espaciados irregularmente y crear su cuadradillo restringido (0.29 segundos de tiempo transcurrido):

Aquíestáelcódigoparaproducirestasparcelas.ExplotaelpolimorfismodeR:sellamaráapoints.quadtreecadavezquelafunciónpointsseapliqueaunobjetoquadtree,porejemplo.Elpoderdeestoesevidenteenlaextremasimplicidaddelafunciónparacolorearlospuntosdeacuerdoconsuidentificadordegrupo:

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

Eltrazadodelacuadrículaensíesunpocomáscomplicadoporquerequiereunrecorterepetidodelosumbralesutilizadosparalaparticióndequadtree,peroelmismoenfoquerecursivoessimpleyelegante.Useunavarianteparaconstruirrepresentacionespoligonalesdeloscuadrantessilodesea.

lines.quadtree<-function(q,xylim,...){i<-q$indexj<-3-q$indexclip<-function(xylim.clip,i,upper){if(upper)xylim.clip[1,i]<-max(q$threshold,xylim.clip[1,i])elsexylim.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,...){}#Nothingtodoatleaves!

Comootroejemplo,generé1,000,000puntosylosdividíengruposde5-9cadauno.Eltiempofuede91.7segundos.

n<-25000#Pointsperclustern.centers<-40#Numberofclustercenterssd<-1/2#Standarddeviationofeachclusterset.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<-5system.time(qt<-quadtree(xy,k))##Setuptomapthefullextentofthequadtree.#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=".")

ComoejemplodecómointeractuarconunSIG,escribamostodaslasceldasdequadtreecomounarchivodeformadepolígonousandolabibliotecashapefiles.Elcódigoemulalasrutinasderecortedelines.quadtree,peroestaveztienequegenerardescripcionesdevectoresdelasceldas.Estossegenerancomomarcosdedatosparausarconlabibliotecashapefiles.

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 usando read.shp o 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 la extensión deseada para xylim aquí para abrir una ventana en una subregión o para expandir la asignación a una región más grande; este código se establece de manera predeterminada en la medida de los puntos).

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

    
respondido por el whuber 21.08.2012 - 17:45
6

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

  1. comienza con una cuadrícula regular
  2. si el polígono tiene menos del umbral, fusione con el vecino alternando (E, S, W, N) en espiral hacia la derecha.
  3. si el polígono tiene menos del umbral, vaya a 2, o vaya al siguiente polígono

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

    
respondido por el Paulo Scardine 21.08.2012 - 03:51
5

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

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

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

    
respondido por el MerseyViking 21.08.2012 - 10:23
1

Otro enfoque es crear una cuadrícula muy fina y usar el algoritmo max-p. enlace

    
respondido por el fgregg 23.06.2014 - 20:51

Lea otras preguntas en las etiquetas