Wie kann ich ein unregelmäßiges Gitter mit mindestens n Punkten erzeugen?

Ist es bei einer großen Stichprobe (~ 1 Million) ungleichmäßig verteilter Punkte möglich, ein unregelmäßiges Raster (in der Größe, aber möglicherweise auch unregelmäßig in der Form?) Zu generieren, das eine festgelegte Mindestanzahl von n Punkten enthält?

Es ist für mich weniger wichtig, wenn generierte "Zellen" eines solchen Gitters genau n Punkte oder mindestens n Punkte enthalten.

Ich kenne Lösungen wie genvecgrid in ArcGIS oder Create Grid Layer in QGIS / mmgis. Sie alle erstellen jedoch reguläre Gitter, die zu einer Ausgabe mit leeren Zellen (kleineres Problem - ich könnte sie einfach verwerfen) oder Zellen mit Punkteanzahl führen kleiner als n (größeres Problem, da ich eine Lösung zum Zusammenfassen dieser Zellen benötigen würde, wahrscheinlich mit einigen Tools von hier ?).

Ich habe vergeblich gegoogelt und bin offen für kommerzielle (ArcGIS & Extensions) oder kostenlose (Python, PostGIS, R) Lösungen.

Ich sehe, dass MerseyViking einen Quadtree empfohlen hat . Ich wollte dasselbe vorschlagen und um es zu erklären, hier ist der Code und ein Beispiel. Der Code ist in geschrieben R, sollte sich aber leicht beispielsweise nach Python portieren lassen.

Die Idee ist bemerkenswert einfach: Teilen Sie die Punkte ungefähr in der Hälfte in der x-Richtung, dann teilen Sie die beiden Hälften rekursiv entlang der y-Richtung, wobei sich die Richtungen auf jeder Ebene abwechseln, bis keine Aufteilung mehr erwünscht ist.

Da die Absicht darin besteht, tatsächliche Punktpositionen zu verschleiern, ist es nützlich, den Teilungen eine gewisse Zufälligkeit zu verleihen . Ein schneller und einfacher Weg, dies zu tun, besteht darin, an einem Quantil einen kleinen Zufallsbetrag abzuspalten, der von 50% abweicht. Auf diese Weise ist es sehr unwahrscheinlich, dass (a) die Aufteilungswerte mit den Datenkoordinaten übereinstimmen, so dass die Punkte eindeutig in durch die Partitionierung erzeugte Quadranten fallen und (b) Punktkoordinaten nicht präzise aus dem Quadtree rekonstruiert werden können.

Da die Absicht besteht, eine Mindestanzahl kvon Knoten in jedem Quadtree-Blatt beizubehalten , implementieren wir eine eingeschränkte Form von Quadtree. Es werden (1) Clustering-Punkte in Gruppen mit jeweils zwischen kund 2 * k-1 Elementen und (2) Mapping der Quadranten unterstützt.

Dieser RCode erstellt einen Baum von Knoten und Terminalblättern, die nach Klassen unterschieden werden. Die Klassenkennzeichnung beschleunigt die Nachbearbeitung, z. B. das unten gezeigte Plotten. Der Code verwendet numerische Werte für die IDs. Dies funktioniert bis zu einer Tiefe von 52 im Baum (unter Verwendung von Doubles; wenn vorzeichenlose lange Ganzzahlen verwendet werden, beträgt die maximale Tiefe 32). Für tiefere Bäume (die in jeder Anwendung höchst unwahrscheinlich sind, da mindestens k* 2 ^ 52 Punkte beteiligt wären) müssten IDs Zeichenfolgen sein.

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

Es ist zu beachten, dass das rekursive Divide-and-Conquer-Design dieses Algorithmus (und folglich der meisten Nachverarbeitungsalgorithmen) bedeutet, dass der Zeitbedarf O (m) und die RAM-Auslastung O (n) mist, wobei die Anzahl von ist Zellen und nist die Anzahl der Punkte. mist proportional zu ndividiert durch die minimalen Punkte pro Zelle,k. Dies ist nützlich, um die Berechnungszeiten abzuschätzen. Wenn zum Beispiel 13 Sekunden benötigt werden, um n = 10 ^ 6 Punkte in Zellen mit 50-99 Punkten (k = 50) zu unterteilen, ist m = 10 ^ 6/50 = 20000. Wenn Sie stattdessen bis 5-9 unterteilen möchten Punkte pro Zelle (k = 5), m ist 10-mal größer, so dass das Timing auf ca. 130 Sekunden ansteigt. (Da der Vorgang des Aufteilens eines Satzes von Koordinaten um ihre Mitten schneller wird, wenn die Zellen kleiner werden, betrug das tatsächliche Timing nur 90 Sekunden.) Bis zu k = 1 Punkt pro Zelle dauert es ungefähr sechsmal länger Noch, oder neun Minuten, und wir können davon ausgehen, dass der Code tatsächlich etwas schneller ist.

Bevor wir fortfahren, generieren wir einige interessante Daten mit unregelmäßigen Abständen und erstellen ihren eingeschränkten Quadtree (0,29 Sekunden verstrichene Zeit):

Hier ist der Code, um diese Diagramme zu erstellen. Der RPolymorphismus points.quadtreewird ausgenutzt : Wird aufgerufen, wenn die pointsFunktion beispielsweise auf ein quadtreeObjekt angewendet wird . Die Stärke davon zeigt sich in der extrem einfachen Funktion, die Punkte gemäß ihrer Clusterkennung einzufärben:

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

Das Zeichnen des Rasters selbst ist etwas komplizierter, da die für die Quadtree-Partitionierung verwendeten Schwellenwerte wiederholt abgeschnitten werden müssen. Derselbe rekursive Ansatz ist jedoch einfach und elegant. Verwenden Sie bei Bedarf eine Variante, um polygonale Darstellungen der Quadranten zu erstellen.

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!

Als weiteres Beispiel habe ich 1.000.000 Punkte generiert und diese in Gruppen von jeweils 5-9 aufgeteilt. Das Timing betrug 91,7 Sekunden.

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

Als Beispiel für die Interaktion mit einem GIS schreiben wir alle Quadtree-Zellen mithilfe der shapefilesBibliothek als Polygon-Shapefile aus . Der Code emuliert die Clipping-Routinen von lines.quadtree, muss jedoch dieses Mal Vektorbeschreibungen der Zellen generieren. Diese werden als Datenrahmen zur Verwendung mit der shapefilesBibliothek ausgegeben .

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

Die Punkte selbst können direkt mithilfe read.shpoder durch Importieren einer Datendatei mit (x, y) Koordinaten gelesen werden .

Anwendungsbeispiel:

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)

(Verwenden Sie eine beliebige Ausdehnung, um xylimhier ein Fenster in eine Teilregion zu öffnen oder die Zuordnung auf eine größere Region zu erweitern. Der Standardwert für diesen Code ist die Ausdehnung der Punkte.)

Dies allein reicht aus: Eine räumliche Verknüpfung dieser Polygone mit den ursprünglichen Punkten identifiziert die Cluster. Einmal identifiziert, erzeugen Datenbank "Zusammenfassungs" -Operationen eine Zusammenfassungsstatistik der Punkte in jeder Zelle.

Überprüfen Sie, ob dieser Algorithmus genügend Anonymität für Ihre Datenprobe bietet:

1. Beginnen Sie mit einem regelmäßigen Raster
2. Wenn das Polygon kleiner als der Schwellenwert ist, verschmelzen Sie mit dem Nachbarn abwechselnd (E, S, W, N) und drehen Sie die Spirale im Uhrzeigersinn.
3. Wenn das Polygon kleiner als der Schwellenwert ist, fahren Sie mit 2 fort, andernfalls fahren Sie mit dem nächsten Polygon fort

Wenn zum Beispiel der Mindestschwellenwert 3 ist:

Wie wäre es, ähnlich wie bei Paulos interessanter Lösung, mit einem Quad-Tree-Unterteilungsalgorithmus?

Stellen Sie die Tiefe ein, in die der Quadtree gehen soll. Sie können auch eine minimale oder maximale Anzahl von Punkten pro Zelle festlegen, sodass einige Knoten tiefer / kleiner als andere sind.

Unterteilen Sie Ihre Welt, indem Sie leere Knoten verwerfen. Spülen und wiederholen, bis die Kriterien erfüllt sind.

Ein anderer Ansatz besteht darin, ein sehr feines Gitter zu erstellen und den max-p-Algorithmus zu verwenden. http://pysal.readthedocs.org/en/v1.7/library/region/maxp.html