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

20

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

radek
źródło
1
Jak „regularna” musi być siatka? Zastanawiam się, czy możesz wykonać jakieś hierarchiczne grupowanie, a następnie po prostu wyciąć dendrogram, aby spełnić twoje potrzeby (chociaż prawdopodobnie rozciąga to, co można by określić jako zwykłą konfigurację przestrzenną). Dokumentacja CrimeStat zawiera kilka dobrych przykładów tego typu grupowania .
Andy W
5
Czy możesz wyjaśnić dokładnie, co rozumiesz przez „nieregularną siatkę”? Brzmi jak oksymoron :-). Co więcej, jaki byłby cel tego ćwiczenia? Zauważ też, że prawdopodobnie potrzebne są dodatkowe kryteria lub ograniczenia: w końcu, jeśli narysujesz kwadrat wokół wszystkich 1 miliona punktów, można go uznać za część siatki i zawierałby więcej niż n z nich. Prawdopodobnie nie zainteresowałbyś się tym trywialnym rozwiązaniem: ale dlaczego nie?
whuber
@AndyW Dzięki. Dobry pomysł i warty odkrycia. Spojrzę. Rozmiar i kształt „siatki” ma drugorzędne znaczenie dla mnie - priorytet (ze względu na poufność danych) jest „hide” n dysponuje za jednego
Radek
@ whuber Również dziękuję. Zgadzam się - ale nie byłem pewien, jak inaczej nazwać takie partycjonowanie. Jak wspomniano powyżej - moją główną motywacją jest prywatność danych. Mając pięć punktów lokalizacji (których nie mogę pokazać na ostatecznej mapie) chciałbym przedstawić je według obszaru, który je obejmuje; i uzyskaj średnią / medianę / etc. wartość za to. Zgadzam się, że byłoby możliwe narysowanie jednego prostokąta lub wypukłego kadłuba reprezentującego je wszystkie - to byłaby chyba najlepsza ochrona prywatności danych? ;] Jednak - bardziej przydatne byłoby przedstawienie go za pomocą granic kształtów, powiedzmy 10 funkcji. Następnie - nadal mogę zachować wzór przestrzenny.
radek,
1
IMO podało twój opis Użyłbym pewnego rodzaju interpolacji i wyświetliłbym mapę rastrową (być może wystarczająca byłaby szerokość pasma wielkości minimalnego N, aby wygładzić dane). Jeśli chodzi o CrimeStat, największe pliki, z których korzystałem, to około 100 000 przypadków, które (moim zdaniem z pewnością zajęłyby trochę czasu). Prawdopodobnie możesz dokonać wstępnej uogólnienia swoich danych, aby przedstawić je jako mniej przypadków i nadal uzyskać pożądane wyniki dla wszystkiego, czego chcesz. To jest naprawdę prosty program, proponuję poświęcić kilka minut, aby go wypróbować i przekonać się sam.
Andy W

Odpowiedzi:

26

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

Quadtree

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

wprowadź opis zdjęcia tutaj


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.

Whuber
źródło
Łał! Fantastyczny. Spróbuję z moimi danymi po powrocie do biura =)
radek
4
Najlepsza odpowiedź @whuber! +1
MerseyViking
1
(1) Możesz czytać pliki kształtów bezpośrednio za pomocą ( między innymi ) shapefilespakietu lub możesz eksportować współrzędne (x, y) w tekście ASCII i czytać je za pomocą read.table. (2) Zalecam zapisanie qtw dwóch formach: po pierwsze, jako punktowy plik kształtu, w xyktórym idpola są zawarte jako identyfikatory klastrów; po drugie, gdzie narysowane segmenty linii lines.quadtreesą zapisywane jako plik kształtu polilinii (lub gdzie analogiczne przetwarzanie zapisuje komórki jako plik kształtu wieloboku). Jest to tak proste, jak modyfikowanie lines.quadtree.leafdanych wyjściowych xylimjako prostokąta. (Zobacz zmiany.)
whuber
1
@ Whubber Bardzo dziękuję za aktualizację. Wszystko działało sprawnie. Zasłużony +50, choć teraz myślę, że zasługuje na +500!
radek
1
Podejrzewam, że z jakiegoś powodu obliczone identyfikatory nie były unikalne. Wprowadź następujące zmiany w definicji quad: (1) zainicjuj id=1; (2) zmiana id/2na id*2w lower=linii; (3) wprowadź podobną zmianę jak id*2+1w upper=wierszu. (Przeredaguję moją odpowiedź, aby to odzwierciedlić.) To powinno również zadbać o obliczenie obszaru: w zależności od Twojego GIS wszystkie obszary będą dodatnie lub wszystkie będą ujemne. Jeśli wszystkie są negatywne, odwróć listy dla xi yw cell.quadtree.leaf.
whuber
6

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:

algorytm

Paulo Scardine
źródło
1
Diabeł tkwi w szczegółach: wydaje się, że takie podejście (lub prawie każdy algorytm skupiania aglomeracyjnego) grozi pozostawieniem małych „osieroconych” punktów rozrzuconych po całym miejscu, które następnie nie mogą zostać przetworzone. Nie twierdzę, że takie podejście jest niemożliwe, ale utrzymałbym zdrowy sceptycyzm przy braku rzeczywistego algorytmu i przykładu jego zastosowania do realistycznego zbioru danych punktowych.
whuber
Rzeczywiście takie podejście może być problematyczne. Jedno zastosowanie tej metody, o którym myślałem, wykorzystuje punkty jako reprezentacje budynków mieszkalnych. Myślę, że ta metoda sprawdziłaby się dobrze w gęsto zaludnionych obszarach. Jednak wciąż zdarzałyby się przypadki, gdy dosłownie jeden lub dwa budynki znajdują się „w szczerym polu” i wymagałoby to wielu iteracji i spowodowałoby, że naprawdę duże obszary ostatecznie osiągnęłyby minimalny próg.
radek
5

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.

MerseyViking
źródło
Dzięki. Jakie oprogramowanie byś polecił?
radek
1
Zasadniczo jest to dobry pomysł. Ale jak powstałyby puste węzły, gdybyś nigdy nie dopuścił mniej niż dodatnią minimalną liczbę punktów na komórkę? (Istnieje wiele rodzajów czworokątów, więc możliwość pustych węzłów wskazuje, że masz na myśli taki, który nie jest dostosowany do danych, co budzi obawy o jego przydatność do zamierzonego zadania.)
whuber
1
Wyobrażam to w ten sposób: wyobraź sobie, że węzeł ma więcej niż maksymalny próg punktów, ale są one skupione w lewym górnym rogu węzła. Węzeł będzie podzielony, ale prawy dolny węzeł potomny będzie pusty, więc można go przyciąć.
MerseyViking
1
Widzę, co robisz (+1). Sztuczka polega na podziale w punkcie określonym przez współrzędne (takie jak ich mediana), co gwarantuje brak pustych komórek. W przeciwnym razie czworokąt określa przede wszystkim przestrzeń zajmowaną przez punkty, a nie same punkty; Twoje podejście staje się skutecznym sposobem realizacji ogólnej idei zaproponowanej przez @Paulo.
whuber