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

20

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.

radek
quelle
1
Wie "normal" muss das Raster sein? Ich frage mich, ob Sie hierarchische Cluster erstellen und dann das Dendrogramm entsprechend Ihren Anforderungen zuschneiden können (obwohl dies wahrscheinlich die Definition einer regulären räumlichen Konfiguration übersteigt). Die CrimeStat-Dokumentation enthält einige gute Beispiele für diese Art von Clustering .
Andy W
5
Könnten Sie bitte genau erklären, was Sie unter einem "unregelmäßigen Gitter" verstehen? Das klingt nach einem Oxymoron :-). Genauer gesagt, was wäre der Zweck dieser Übung? Beachten Sie auch, dass wahrscheinlich zusätzliche Kriterien oder Einschränkungen erforderlich sind: Wenn Sie ein Quadrat um alle 1 Million Punkte zeichnen, könnte es als Teil eines Gitters betrachtet werden und mehr als n davon enthalten. Sie würden sich wahrscheinlich nicht für diese triviale Lösung interessieren: aber warum nicht genau?
whuber
@ AndyW Danke. Gute Idee und es lohnt sich, sie zu erkunden. Werde mal schauen. Größe & Form von ‚Grid‘ ist von untergeordneter Bedeutung für mich - Priorität (aufgrund des Datenschutzes) ist zu ‚verstecken‘ n kennzeichnet hinter einem
radek
@whuber Danke auch. Ich stimme dem zu - war mir aber nicht sicher, wie ich eine solche Unterteilung sonst nennen könnte. Wie oben erwähnt - meine Hauptmotivation ist der Datenschutz. Mit fünf Punktpositionen (die ich auf der endgültigen Karte nicht anzeigen kann) möchte ich sie nach bedeckendem Gebiet darstellen. und erhalte mean / median / etc. Wert dafür. Ich bin damit einverstanden, dass es möglich wäre, ein Rechteck oder eine konvexe Hülle zu zeichnen, die alle darstellen - das wäre der ultimative Datenschutz, denke ich. ;] Es wäre jedoch sinnvoller, es durch Begrenzung von Formen darzustellen, sagen wir 10 Merkmale. Dann kann ich noch räumliche Muster bewahren.
Radek
1
Wenn Sie mir eine Beschreibung geben, würde ich eine Art Interpolation verwenden und eine Rasterkarte anzeigen (möglicherweise würde eine adaptive Bandbreite von der Größe Ihres minimalen N ausreichen, um die Daten zu glätten). Was CrimeStat anbelangt, so waren die größten Dateien, die ich verwendet habe, meiner Meinung nach etwa 100.000 Fälle (und das Clustering würde sicherlich einige Zeit in Anspruch nehmen). Es ist wahrscheinlich, dass Sie Ihre Daten vorab generalisieren können, um sie in weniger Fällen darzustellen, und dennoch wünschenswerte Ergebnisse für alles zu erzielen, was Sie möchten. Es ist ein wirklich einfaches Programm, ich würde vorschlagen, sich nur ein paar Minuten Zeit zu nehmen, um es auszuprobieren und sich selbst davon zu überzeugen.
Andy W

Antworten:

26

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

Quadtree

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

Bildbeschreibung hier eingeben


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.

whuber
quelle
Wow! Fantastisch. Ich werde es mit meinen Daten
versuchen,
4
Top Antwort @whuber! +1
MerseyViking
1
(1) Sie können Shapefiles direkt mit ( unter anderem ) dem shapefilesPaket lesen oder Sie können (x, y) Koordinaten in ASCII-Text exportieren und mit lesen read.table. (2) Ich empfehle, qtin zwei Formen zu schreiben : erstens als Punkt-Shapefile, xyin dem die idFelder als Cluster-IDs enthalten sind; Zweitens, wenn die durchgezeichneten Liniensegmente lines.quadtreeals Polylinien-Shapefile ausgeschrieben werden (oder wenn bei analoger Verarbeitung die Zellen als Polygon-Shapefile geschrieben werden). Dies ist so einfach wie das Ändern lines.quadtree.leafder Ausgabe xylimals Rechteck. (Siehe die Änderungen.)
whuber
1
@whubber Vielen Dank für ein Update. Alles hat reibungslos funktioniert. Gut verdient +50, obwohl ich jetzt denke, dass es +500 verdient!
Radek
1
Ich vermute, die berechneten IDs waren aus irgendeinem Grund nicht eindeutig. Nehmen Sie diese Änderungen in der Definition vor quad: (1) initialisieren id=1; (2) Änderung id/2zu id*2der lower=Leitung; (3) Nehmen Sie eine ähnliche Änderung wie id*2+1in der upper=Zeile vor. (Ich werde meine Antwort entsprechend anpassen.) Dies sollte auch die Flächenberechnung berücksichtigen: Abhängig von Ihrem GIS sind alle Flächen positiv oder alle negativ. Wenn sie alle negativ sind, kehren Sie die Listen für xund yin um cell.quadtree.leaf.
Whuber
6

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

Algorithmus

Paulo Scardine
quelle
1
Der Teufel steckt im Detail: Dieser Ansatz (oder fast jeder agglomerative Clustering-Algorithmus) droht, überall kleine "Orphan" -Punkte zu hinterlassen, die dann nicht verarbeitet werden können. Ich sage nicht, dass dieser Ansatz unmöglich ist, aber ich würde ohne einen tatsächlichen Algorithmus und ein Beispiel für dessen Anwendung auf einen realistischen Punktdatensatz eine gesunde Skepsis aufrechterhalten.
Whuber
In der Tat könnte dieser Ansatz problematisch sein. Eine Anwendung dieser Methode, über die ich nachgedacht habe, verwendet Punkte als Darstellungen von Wohngebäuden. Ich denke, diese Methode würde in dichter besiedelten Gebieten gut funktionieren. Es würde jedoch immer noch Fälle geben, in denen sich buchstäblich ein oder zwei Gebäude „mitten im Nirgendwo“ befinden und viele Iterationen erforderlich sind. Dies würde dazu führen, dass wirklich große Flächen die Mindestschwelle erreichen.
Radek
5

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.

MerseyViking
quelle
Vielen Dank. Welche Software würdest du dafür empfehlen?
Radek
1
Im Prinzip ist das eine gute Idee. Aber wie würden leere Knoten entstehen, wenn Sie niemals weniger als eine positive Mindestanzahl von Punkten pro Zelle zulassen? (Da es viele Arten von Quadtrees gibt, weist die Möglichkeit leerer Knoten darauf hin, dass Sie eines im Sinn haben, das nicht an die Daten angepasst ist, was Bedenken hinsichtlich seiner Nützlichkeit für die beabsichtigte Aufgabe aufwirft.)
whuber
1
Ich stelle es mir so vor: Stellen Sie sich vor, ein Knoten hat mehr als die maximale Punkteschwelle, aber sie sind oben links im Knoten gruppiert. Der Knoten wird unterteilt, der untergeordnete Knoten unten rechts ist jedoch leer, sodass er beschnitten werden kann.
MerseyViking
1
Ich sehe, was du tust (+1). Der Trick besteht darin, an einem Punkt zu unterteilen, der durch die Koordinaten bestimmt wird (wie z. B. deren Median), wodurch garantiert wird, dass keine leeren Zellen vorhanden sind. Ansonsten wird der Quadtree in erster Linie durch den Platz bestimmt, den die Punkte einnehmen und nicht die Punkte selbst; Ihr Ansatz wird dann zu einer effektiven Möglichkeit, die von @Paulo vorgeschlagene generische Idee umzusetzen.
whuber