Veo MerseyViking ha recomendado un quadtree. Yo iba a sugerir lo mismo y con el fin de explicar, aquí está el código y un ejemplo. El código está escrito en R pero debería puerto fácilmente, digamos, Python.
La idea es muy sencilla: dividir los puntos aproximadamente en la mitad en la dirección x, entonces recursivamente divide las dos mitades a lo largo de la dirección del eje y, alternando las direcciones en cada nivel, hasta que no hay más división que se desea.
Debido a que la intención es la de disimular la real ubicaciones de punto, es útil introducir un poco de aleatoriedad en las divisiones. Una rápida manera simple de hacer esto es para dividir a un cuantil establecer una pequeña cantidad aleatoria lejos del 50%. En esto de la moda (a) la división de valores es muy poco probable que coinciden con las coordenadas de los datos, de modo que los puntos se caen de forma exclusiva en cuadrantes creado por la partición, y (b) las coordenadas del punto será imposible reconstruir con precisión desde el quadtree.
Porque la intención es mantener una mínima cantidad k de nodos dentro de cada quadtree de la hoja, poner en práctica un restringido de quadtree. Se va a apoyar (1) de la agrupación de los puntos en grupos de entre k y 2*k-1 elementos de cada uno y (2) la asignación de los cuadrantes.
Este R código crea un árbol de nodos y terminales de las hojas, que la distingue por clase. La clase etiquetado acelera de post-procesamiento, tales como trazado, que se muestra a continuación. El código utiliza valores numéricos para los identificadores. Esto funciona hasta profundidades de 52 en el árbol (el uso de dobles; si unsigned long enteros se utilizan, la profundidad máxima es de 32). Para profundizar en los árboles (el cual es altamente improbable en cualquier aplicación, ya que, al menos, k * 2^52 puntos estarían involucrados), ids tendría que 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 recurrente dividir y conquistar el diseño de este algoritmo (y, en consecuencia, de la mayoría de los post-procesamiento de los algoritmos) significa que el requisito de tiempo es O(m) y el uso de la RAM es O(n) donde m es el número de células y n es el número de puntos. m es proporcional a n dividido por el mínimo de puntos por celda, k. Esto es útil para la estimación de los tiempos de cómputo. Por ejemplo, si se tarda 13 segundos para la partición de n=10^6 puntos en las células de 50 a 99 puntos (k=50), m = 10^6/50 = 20000. Si desea lugar a la partición de abajo a 5-9 puntos por celda (k=5), m es 10 veces mayor, por lo que el calendario va hasta alrededor de 130 segundos. (Debido a que el proceso de dividir un conjunto de coordenadas alrededor de sus medios se hace más rápido a medida que las células se hacen más pequeños, el momento en que se fue de sólo 90 segundos). Ir todo el camino a k=1 punto por cada celda, se va a tomar alrededor de seis veces más todavía, o nueve minutos, y se puede esperar que el código, en realidad para ser un poco más rápido que eso.
Antes de continuar, vamos a generar unos interesantes irregularmente espaciados de datos y crear sus restringido quadtree (0.29 segundos de tiempo transcurrido):
![Quadtree]()
Aquí está el código para producir estas parcelas. Explota R's polimorfismo: points.quadtree será llamado cada vez que el points función se aplica a un quadtree objeto, por ejemplo. El poder de esto es evidente en la extrema simplicidad de la función de color de los puntos de acuerdo a su identificador de clúster:
points.quadtree <- function(q, ...) {
points(q$lower, ...); points(q$upper, ...)
}
points.quadtree.leaf <- function(q, ...) {
points(q$value, col=hsv(q$id), ...)
}
El trazado de la red propia es un poco más complicada porque requiere repetidas recorte de los umbrales utilizados para el quadtree de partición, pero el mismo recursiva enfoque es simple y elegante. El uso de una variante para la construcción de planta poligonal, las representaciones de los cuadrantes, si se 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, yo generado 1.000.000 de puntos y dividido en grupos de 5 a 9 de 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=".")
![enter image description here]()
Como un ejemplo de cómo interactuar con un SIG, vamos a escribir todos los quadtree células como un shapefile de polígonos utilizando la shapefiles biblioteca. El código emula las rutinas de recorte de lines.quadtree, pero esta vez se tiene para generar el vector de las descripciones de las células. Estos son de salida como de los marcos de datos para su uso con el shapefiles biblioteca.
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 se pueden leer directamente el uso de read.shp o mediante la importación de un archivo de datos de (x,y) las coordenadas.
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)
(El uso de cualquier medida deseada para xylim aquí a la ventana en una subregión o ampliar la asignación de una región más grande; este código de valores predeterminados para la medida de los puntos).
Esto por sí solo es suficiente: una unión espacial de los polígonos a los puntos originales identificará los grupos. Una vez identificado, la base de datos de "resumir" operaciones de generar las estadísticas de resumen de los puntos dentro de cada célula.
