Jak wygenerować nieregularną siatkę zawierającą minimum n punktów?

Biorąc pod uwagę dużą (~ 1 milion) próbkę nierównomiernie rozmieszczonych punktów - czy można wygenerować nieregularną siatkę (pod względem wielkości, ale może również mieć nieregularny kształt, jeśli to możliwe?), Która będzie zawierać określoną minimalną ilość n punktów?

Dla mnie nie ma większego znaczenia, jeśli generetowane „komórki” takiej siatki zawierają dokładnie n liczby punktów lub co najmniej n punktów.

Jestem świadomy rozwiązań takich jak genvecgrid w ArcGIS lub Create Grid Layer w QGIS / mmgis, jednak wszystkie one utworzą regularne siatki, które skutkują wyjściem z pustymi komórkami (mniejszy problem - mógłbym je po prostu odrzucić) lub komórkami z liczbą punktów mniej niż n (większy problem, ponieważ potrzebuję rozwiązania do agregacji tych komórek, prawdopodobnie przy użyciu niektórych narzędzi stąd ?).

Googling bezskutecznie i jestem otwarty na rozwiązania komercyjne (ArcGIS i rozszerzenia) lub bezpłatne (Python, PostGIS, R).

Widzę MerseyViking zaleciła QuadTree . Chciałem zasugerować to samo i aby to wyjaśnić, oto kod i przykład. Kod jest zapisany, Rale powinien być łatwy do przeniesienia, powiedzmy, do Pythona.

Pomysł jest niezwykle prosty: podziel punkty w przybliżeniu na pół w kierunku x, a następnie rekurencyjnie podziel dwie połówki wzdłuż kierunku y, zmieniając kierunki na każdym poziomie, aż nie będzie już konieczne dzielenie.

Ponieważ celem jest ukrycie rzeczywistych lokalizacji punktów, przydatne jest wprowadzenie losowości do podziałów . Jednym szybkim i prostym sposobem na dokonanie tego jest podzielenie na zestaw kwantyli małej losowej kwoty od 50%. W ten sposób (a) jest mało prawdopodobne, aby wartości podziału były zbieżne ze współrzędnymi danych, tak że punkty będą jednoznacznie wpadać do kwadrantów utworzonych przez podział, oraz (b) współrzędnych punktów nie będzie można dokładnie odtworzyć z kwadratu.

Ponieważ intencją jest utrzymanie minimalnej kliczby węzłów w każdym liściu czterokąta, implementujemy ograniczoną formę czterokąta. Będzie obsługiwał (1) punkty grupowania w grupy posiadające pomiędzy ki 2 * k-1 każdy element oraz (2) mapowanie kwadrantów.

Ten Rkod tworzy drzewo węzłów i liści końcowych, rozróżniając je według klasy. Etykietowanie klas przyspiesza przetwarzanie końcowe, takie jak drukowanie, przedstawione poniżej. Kod używa wartości liczbowych dla identyfikatorów. Działa to do głębokości 52 w drzewie (przy użyciu podwójnych; jeśli używane są długie liczby całkowite bez znaku, maksymalna głębokość wynosi 32). W przypadku głębszych drzew (które są bardzo mało prawdopodobne w żadnej aplikacji, ponieważ w grę wchodzą co najmniej k* 2 ^ 52 punkty), identyfikatory musiałyby być łańcuchami.

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

Należy zauważyć, że rekurencyjna konstrukcja tego algorytmu typu „dziel i rządź” (i w konsekwencji większość algorytmów przetwarzania końcowego) oznacza, że ​​wymagany czas to O (m), a użycie pamięci RAM to O (n), gdzie mjest liczba komórki i njest liczbą punktów. mjest proporcjonalne don podzielonej przez minimalną punktów na komórkę,k. Jest to przydatne do szacowania czasów obliczeń. Na przykład, jeśli podzielenie n = 10 ^ 6 punktów na komórki 50-99 punktów (k = 50) zajmuje 13 sekund, m = 10 ^ 6/50 = 20000. Jeśli zamiast tego chcesz podzielić na 5-9 punktów na komórkę (k = 5), m jest 10 razy większy, więc czas wzrasta do około 130 sekund. (Ponieważ proces dzielenia zestawu współrzędnych wokół ich środków przyspiesza, gdy komórki stają się mniejsze, faktyczny czas wyniósł zaledwie 90 sekund.) Aby przejść do k = 1 punkt na komórkę, zajmie to około sześć razy dłużej jeszcze dziewięć minut i możemy spodziewać się, że kod będzie trochę szybszy.

Zanim przejdziemy dalej, wygenerujmy interesujące dane o nieregularnych odstępach i stwórzmy ograniczony quadree (czas, który upłynął 0,29 sekundy):

Oto kod do tworzenia tych wykresów. Wykorzystuje Rpolimorfizm: points.quadtreebędzie wywoływany za każdym razem, gdy pointsfunkcja zostanie zastosowana na przykład do quadtreeobiektu. Siła tego jest widoczna w ekstremalnej prostocie funkcji pokolorowania punktów zgodnie z ich identyfikatorem skupienia:

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

Wykreślenie samej siatki jest nieco trudniejsze, ponieważ wymaga wielokrotnego przycinania progów użytych do partycjonowania quadtree, ale to samo rekurencyjne podejście jest proste i eleganckie. W razie potrzeby użyj wariantu, aby zbudować wielokątne reprezentacje ćwiartek.

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!

Jako inny przykład wygenerowałem 1 000 000 punktów i podzieliłem je na grupy po 5-9 sztuk. Czas wyniósł 91,7 sekundy.

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

Jako przykład interakcji z GIS , wypiszmy wszystkie komórki z czterema drzewami jako plik kształtu wielokąta za pomocą shapefilesbiblioteki. Kod emuluje procedury obcinania lines.quadtree, ale tym razem musi wygenerować opisy wektorowe komórek. Są one wyprowadzane jako ramki danych do użytku z shapefilesbiblioteką.

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

Same punkty można odczytać bezpośrednio, używając read.shplub importując plik danych o współrzędnych (x, y).

Przykład zastosowania:

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)

(Użyj tutaj dowolnego pożądanego zasięgu, aby xylimprzejść do podregionu lub rozszerzyć mapowanie na większy region; ten kod domyślnie określa zakres punktów).

Już samo to wystarczy: przestrzenne połączenie tych wielokątów z pierwotnymi punktami zidentyfikuje skupiska. Po zidentyfikowaniu operacje „podsumowania” w bazie danych wygenerują statystyki podsumowujące punkty w każdej komórce.

Sprawdź, czy ten algorytm zapewnia wystarczającą anonimowość dla próbki danych:

1. zacznij od zwykłej siatki
2. jeśli wielokąt ma mniej niż próg, połącz z sąsiadującym naprzemiennie (E, S, W, N) spiralnie zgodnie z ruchem wskazówek zegara.
3. jeśli wielokąt ma mniej niż próg, przejdź do 2, w przeciwnym razie przejdź do następnego wielokąta

Na przykład, jeśli minimalny próg wynosi 3:

Podobnie jak w przypadku interesującego rozwiązania Paulo, co powiesz na zastosowanie algorytmu podziału drzewa quad?

Ustaw głębokość, na którą ma iść quadtree. Możesz również mieć minimalną lub maksymalną liczbę punktów na komórkę, aby niektóre węzły były głębsze / mniejsze niż inne.

Podziel swój świat, odrzucając puste węzły. Opłucz i powtarzaj, aż kryteria zostaną spełnione.

Innym podejściem jest utworzenie bardzo drobnej siatki i użycie algorytmu max-p. http://pysal.readthedocs.org/en/v1.7/library/region/maxp.html