Określanie, czy punkt jest otoczony za pomocą przetwarzania rastrowego

Usiłuję ulepszyć niezwykle uciążliwy proces wektor / python dla modelu naturalnego zagrożenia. W tej chwili mamy długi skrypt, który generuje linie odległości / namiaru z danego punktu, aby określić:

1. rodzaj wielokąta, który przecina (np. las, trawa, bagna itp.)
2. odległość do tego wielokąta
3. ile z tych linii przecina wielokąty, aby określić, jak „otoczone” jest.

Jest o wiele więcej zaangażowanych, ale to jest sedno. Usiłuję znaleźć sposób, aby to poprawić i obecnie jestem zaskoczony częścią 3. Pomysł polega na ustaleniu, czy punkt jest całkowicie otoczony wielokątami, na przykład w odległości 200 m

Tak więc na załączonym obrazie chciałbym, aby punkt A był oznaczony jako bardziej zagrożony niż punkt B, ponieważ jest on całkowicie otoczony przez moje wielokąty. Powtarza się to dla około 13 milionów punktów, więc nie jest to małe zadanie i wolę mieć powierzchnię, z której można czerpać wartości, niż uruchamiać nasz skrypt. Myślę, że musi istnieć różnorodność narzędzi hydrologicznych lub ścieżka kosztów, aby to zrobić, ale wydaje mi się, że nie mogę się tym zająć.

Jak mogłem o tym poradzić?

Istnieje rozwiązanie ścieżki kosztów, ale będziesz musiał go sam kodować. Oto, jak mogłoby to wyglądać po zastosowaniu do każdego punktu na obrazie w pytaniu (nieco zgrubnym, aby przyspieszyć obliczenia):

Czarne komórki są częściami otaczających wielokątów. Kolory, od jasnopomarańczowego (krótkiego) do niebieskiego (długiego), pokazują maksymalną odległość (maksymalnie do 50 komórek), którą można osiągnąć poprzez przejście przez linię widzenia bez przechwytywania komórek wielokąta. (Każda komórka poza zakresem tego obrazu jest traktowana jako część wielokątów).

Omówmy skuteczny sposób, aby to zrobić za pomocą rastrowej reprezentacji danych. W tej reprezentacji wszystkie „otaczające” wielokątne komórki będą miały, powiedzmy, niezerowe wartości, a każda komórka, którą można „przejrzeć”, będzie miała wartość zerową.

Krok 1: Wstępne obliczenie struktury danych sąsiedztwa

Najpierw musisz zdecydować, co to znaczy blokować jedną komórkę. Jedną z najpiękniejszych zasad, jakie mogę znaleźć, jest taka: używając współrzędnych całkowitych dla wierszy i kolumn (i zakładając komórki kwadratowe), zastanówmy się, które komórki mogą blokować komórkę (i, j) z widoku na początku (0,0). Nominuję komórkę (i ', j'), która jest najbliżej odcinka linii łączącego (i, j) z (0,0) spośród wszystkich komórek, których współrzędne różnią się od i i co najwyżej o 1. Ponieważ to nie zawsze dać unikalne rozwiązanie (na przykład z (i, j) = (1,2) oba (0,1) i (1,1) będą działać równie dobrze), potrzebne są pewne środki do rozwiązania powiązań. Byłoby miło, gdyby w tej rozdzielczości powiązań szanowano symetrie okrągłych sąsiedztw w siatkach: negowanie albo współrzędnych, albo zmiana współrzędnych zachowuje te sąsiedztwa. Dlatego możemy zdecydować, które komórki blokują (i,

Ilustracją tej reguły jest napisany w poniższym kodzie prototypowym R. Ten kod zwraca strukturę danych, która będzie wygodna do określania „otaczania” dowolnych komórek w siatce.

screen <- function(k=1) {
  #
  # Returns a data structure:
  #   $offset is an array of offsets
  #   $screened is a parallel array of screened offset indexes.
  #   $distance is a parallel array of distances.
  # The first index always corresponds to (0,0).
  #
  screened.by <- function(xy) {
    uv <- abs(xy)
    if (reversed <- uv[2] > uv[1]) {
      uv <- rev(uv)
    }
    i <- which.min(c(uv[1], abs(uv[1]-uv[2]), uv[2]))
    ij <- uv + c(floor((1-i)/3), floor(i/3)-1)
    if (reversed) ij <- rev(ij)
    return(ij * sign(xy))
  }
  #
  # For each lattice point within the circular neighborhood,
  # find the unique lattice point that screens it from the origin.
  #
  xy <- subset(expand.grid(x=(-k:k), y=(-k:k)), 
               subset=(x^2+y^2 <= k^2) & (x != 0 | y != 0))
  g <- t(apply(xy, 1, function(z) c(screened.by(z), z)))
  #
  # Sort by distance from the origin.
  #
  colnames(g) <- c("x", "y", "x.to", "y.to")
  ij <- unique(rbind(g[, 1:2], g[, 3:4]))
  i <- order(abs(ij[,1]), abs(ij[,2])); ij <- ij[i, , drop=FALSE]
  rownames(ij) <- 1:length(i)
  #
  # Invert the "screened by" relation to produce the "screened" relation.
  #
  # (Row, column) offsets.
  ij.df <- data.frame(ij, i=1:length(i))
  #
  # Distances from the origin (in cells).
  distance <- apply(ij, 1, function(u) sqrt(sum(u*u)))
  #
  # "Screens" relation (represented by indexes into ij).
  g <- merge(merge(g, ij.df), ij.df, 
             by.x=c("x.to", "y.to"), by.y=c("x","y"))
  g <- subset(g, select=c(i.x, i.y))
  h <- by(g$i.y, g$i.x, identity)

  return( list(offset=ij, screened=h, distance=distance) )
}

Wartość screen(12)została wykorzystana do stworzenia tego obrazu relacji przesiewania: strzałki wskazują od komórek do tych, które natychmiast je przesiewają. Odcienie są proporcjonalne do odległości do źródła, które znajduje się w środku tego sąsiedztwa:

Obliczenia te są szybkie i należy je wykonać tylko raz dla danej okolicy. Na przykład, gdy patrzysz 200 m na siatce z 5 m komórkami, rozmiar sąsiedztwa wyniesie 200/5 = 40 jednostek.

Krok 2: Zastosowanie obliczeń do wybranych punktów

Reszta jest prosta: aby ustalić, czy komórka znajdująca się w (x, y) (we współrzędnych wiersza i kolumny) jest „otoczona” w odniesieniu do tej struktury danych sąsiedztwa, wykonaj test rekurencyjnie, zaczynając od przesunięcia (i, j) = (0,0) (początek sąsiedztwa). Jeśli wartość w siatce wielokątów w (x, y) + (i, j) jest różna od zera, widoczność jest tam zablokowana. W przeciwnym razie będziemy musieli rozważyć wszystkie przesunięcia, które mogły zostać zablokowane przy przesunięciu (i, j) (które znajdują się w czasie O (1) przy użyciu zwróconej przez dane struktury screen). Jeśli nie ma żadnych zablokowanych, osiągnęliśmy obwód i stwierdziliśmy, że (x, y) nie jest otoczony, więc zatrzymujemy obliczenia (i nie zawracamy sobie głowy sprawdzaniem pozostałych punktów w okolicy).

Możemy gromadzić jeszcze bardziej przydatne informacje, śledząc najdalszą odległość pola widzenia osiągniętą podczas algorytmu. Jeśli jest to mniej niż pożądany promień, komórka jest otoczona; inaczej nie jest.

Oto Rprototyp tego algorytmu. Jest dłuższy, niż się wydaje, ponieważ Rnatywnie nie obsługuje (prostej) struktury stosu wymaganej do zaimplementowania rekurencji, więc stos również musi zostać zakodowany. Właściwy algorytm zaczyna się około dwóch trzecich drogi i potrzebuje tylko kilkunastu linii. (A połowa z nich po prostu radzi sobie z sytuacją wokół krawędzi siatki, sprawdzając indeksy poza zakresem w sąsiedztwie. Można to usprawnić, rozszerzając siatkę wielokątów o krzędy i kolumny wokół jej obwodu, eliminując wszelkie potrzeba sprawdzenia zakresu indeksu kosztem nieco więcej pamięci RAM, aby utrzymać siatkę wielokątów.)

#
# Test a grid point `ij` for a line-of-sight connection to the perimeter
# of a circular neighborhood.  
#   `xy` is the grid.
#   `counting` determines whether to return max distance or count of stack ops.
#   `perimeter` is the assumed values beyond the extent of `xy`.
#
# Grid values of zero admit light; all others block visibility
# Returns maximum line-of-sight distance found within `nbr`.
#
panvisibility <- function(ij, xy, nbr=screen(), counting=FALSE, perimeter=1) {
  #
  # Implement a stack for the algorithm.
  #
  count <- 0 # Stack count
  stack <- list(ptr=0, s=rep(NA, dim(nbr$offset)[1]))
  push <- function(x) {
    n <- length(x)
    count <<- count+n         # For timing
    stack$s[1:n + stack$ptr] <<- x
    stack$ptr <<- stack$ptr+n
  }
  pop <- function() {
    count <<- count+1         # For timing
    if (stack$ptr <= 0) return(NULL)
    y <- stack$s[stack$ptr]
    #stack$s[stack$ptr] <<- NA # For debugging
    stack$ptr <<- stack$ptr - 1
    return(y)
  }
  #
  # Initialization.
  #
  m <- dim(xy)[1]; n <- dim(xy)[2]
  push(1) # Stack the *indexes* of nbr$offset and nbr$screened.
  dist.max <- -1
  #
  # The algorithm.
  #
  while (!is.null(i <- pop())) {
    cell <- nbr$offset[i, ] + ij
    if (cell[1] <= 0 || cell[1] > m || cell[2] <= 0 || cell[2] > n) {
      value <- perimeter
    } else {  
      value <- xy[cell[1], cell[2]]
    }
    if (value==0) {
      if (nbr$distance[i] > dist.max) dist.max <- nbr$distance[i]
      s <- nbr$screened[[paste(i)]]
      if (is.null(s)) {
        #exited = TRUE
        break
      }
      push(s)
    }
  }
  if (counting) return ( count )
  return(dist.max)
}

W tym przykładzie komórki wielokąta są czarne. Kolory podają maksymalną odległość linii widzenia (do 50 komórek) dla komórek niepoligonalnych, od jasnopomarańczowego dla krótkich odległości do ciemnoniebieskiego dla najdłuższych odległości. (Komórki mają szerokość jednej jednostki i wysokość.) Widoczne smugi są tworzone przez małe „wyspy” wielokątów pośrodku „rzeki”: każda blokuje długą linię innych komórek.

Analiza algorytmu

Struktura stosu dokonuje pierwszego wyszukiwania głębokości wykresu widoczności sąsiedztwa w celu znalezienia dowodów, że komórka nie jest otoczona. Jeśli komórki są daleko od dowolnego wielokąta, to wyszukiwanie będzie wymagało sprawdzenia tylko komórek O (k) pod kątem okrągłego sąsiedztwa o promieniu-k. Najgorsze przypadki występują, gdy w sąsiedztwie znajduje się niewielka liczba rozproszonych komórek wielokąta, ale mimo to nie można do końca dotrzeć do granicy sąsiedztwa: wymagają one sprawdzenia prawie wszystkich komórek w każdym sąsiedztwie, czyli O (k ^ 2) operacja.

Następujące zachowanie jest typowe dla tego, co zostanie napotkane. W przypadku małych wartości k, chyba że wielokąty wypełnią większość siatki, większość niepoligonalnych komórek będzie oczywiście nieuziemiona, a algorytm skaluje się jak O (k). W przypadku wartości pośrednich skalowanie zaczyna wyglądać jak O (k ^ 2). Ponieważ k staje się naprawdę duże, większość komórek zostanie otoczona, a fakt ten można ustalić na długo przed sprawdzeniem całego sąsiedztwa: wysiłek obliczeniowy algorytmu osiąga w ten sposób praktyczny limit. Limit ten jest osiągany, gdy promień sąsiedztwa zbliża się do średnicy największych połączonych niepoligonalnych regionów w siatce.

Jako przykład używam countingopcji zakodowanej w prototypie, screenaby zwrócić liczbę operacji stosu używanych w każdym wywołaniu. Mierzy to wysiłek obliczeniowy. Poniższy wykres przedstawia średnią liczbę operacji stosu w funkcji promienia sąsiedztwa. Wykazuje przewidywane zachowanie.

Możemy to wykorzystać do oszacowania obliczeń potrzebnych do oceny 13 milionów punktów na siatce. Załóżmy, że zastosowano sąsiedztwo k = 200/5 = 40. Wtedy średnio potrzeba będzie kilkuset operacji na stosie (w zależności od złożoności siatki wieloboków i miejsca, w którym 13 milionów punktów znajduje się względem wieloboków), co oznacza, że ​​w efektywnie skompilowanym języku, co najwyżej kilka tysięcy prostych operacji numerycznych będą wymagane (dodawanie, mnożenie, odczyt, zapis, przesunięcie itp.). Większość komputerów będzie w stanie ocenić otoczenie około miliona punktów w tym tempie. (TheRimplementacja jest znacznie wolniejsza, ponieważ jest słaba w tego rodzaju algorytmie, dlatego można ją uważać tylko za prototyp). W związku z tym możemy mieć nadzieję, że wydajna implementacja w rozsądnie wydajnym i odpowiednim języku - C ++ i przychodzi mi na myśl Python - może zakończyć ocenę 13 milionów punktów w ciągu minuty lub mniej, zakładając , że cała siatka wielokątów znajduje się w pamięci RAM.

Gdy siatka jest zbyt duża, aby zmieścić się w pamięci RAM, tę procedurę można zastosować do sąsiadujących części siatki. Muszą nakładać się tylko na krzędy i kolumny; podczas mozaikowania wyników weź maksima na zakładki.

Inne aplikacje

„Fetch” z wód jest ściśle związana z „surroundedness” swoich punktów. W rzeczywistości, jeśli użyjemy promienia sąsiedztwa równego lub większego niż średnica akwenu, stworzymy siatkę pobierania (bezkierunkowego) w każdym punkcie akwenu. Dzięki zastosowaniu mniejszego promienia sąsiedztwa uzyskamy przynajmniej dolną granicę pobierania we wszystkich punktach pobierania o najwyższej wartości, co w niektórych aplikacjach może być wystarczające (i może znacznie zmniejszyć wysiłek obliczeniowy). Wariant tego algorytmu, który ogranicza relację „screened by” do określonych kierunków, byłby jednym ze sposobów skutecznego obliczenia pobierania w tych kierunkach. Zauważ, że takie warianty wymagają modyfikacji kodu dla screen; kod dla w panvisibilityogóle się nie zmienia.

Zdecydowanie mogę zobaczyć, jak można to zrobić za pomocą rozwiązania rastrowego, ale biorąc pod uwagę nawet zmniejszoną liczbę punktów, oczekiwałbym bardzo dużej / wysokiej rozdzielczości, a zatem trudnej do przetworzenia siatki lub zestawu siatek. Biorąc to pod uwagę, zastanawiam się, czy wykorzystanie topologii w gdb może być bardziej wydajne. Możesz znaleźć wszystkie puste przestrzenie za pomocą czegoś takiego jak:

arcpy.env.workspace = 'myGDB'
arcpy.CreateTopology_management('myGDB', 'myTopology', '')    
arcpy.AddFeatureClassToTopology_management('myTopology', 'myFeatures', '1','1')    
arcpy.AddRuleToTopology_management ('myToplogy', 'Must Not Have Gaps (Area)', 'myFeatures', '', '', '')    
arcpy.ValidateTopology_management('myTopology', 'Full_Extent')
arcpy.ExportTopologyErrors_management('myTopology', 'myGDB', 'topoErrors')
arcpy.FeatureToPolygon_management('topoErrors_line','topoErrorsVoidPolys', '0.1')`

możesz wtedy pracować z topoErrorsVoidPolysnormalnym wzorcem Intersect_analysis()lub czymkolwiek. Być może będziesz musiał zadzierać z wydobywaniem interesujących cię biegunów topoErrorsVoidPolys. @whuber ma wiele świetnych postów na ten temat w innych miejscach na gis.stackexchange.com.

Jeśli naprawdę chcesz przejść do rastra ... Zrobiłbym coś zgodnie z tym pseudo kodem (nie żałuj tylko dlatego, że to oczywiste, że jestem AML!: P)

1. rasteryzuj punkty („pts_g”) i polys („polys_g” (
2. pustki = regiongroup (con (isnull (polys_g), 1))
3. być może trzeba będzie zrobić coś, aby ulepszyć puste przestrzenie, aby wyeliminować niechciany zewnętrzny wielokąt / otwarty obszar wszechświata
4. pts_surrounded = con (pustki, pts_g)

Po prostu wymyślenie tego, więc może być konieczne udoskonalenie.

Jeśli używasz progowej wartości odległości (tutaj mówisz o 200 m), najlepszym rozwiązaniem jest zastosowanie analizy wektorowej:

1) utwórz bufor o długości 200 m wokół każdego punktu (czarny na ilustracji)

2) użyj narzędzia przecięcia (analizy) między buforem a wielokątami (na niebiesko na ilustracji). Będzie ładniej wyglądać, jeśli zrobisz to między granicami otaczających wielokątów i bufora, ale końcowy wynik jest taki sam.

3) użyj funkcji do wielokąta (zarządzania), aby utworzyć wielokąty, w których twoje punkty są całkowicie otoczone (czerwony na ilustracji)

4) wybierz warstwy według lokalizacji (zarządzanie) lub łączenia przestrzennego (analiza), aby zidentyfikować otoczone punkty. Zastosowanie łączenia przestrzennego pozwala uzyskać informacje o osadzonym wielokącie (obszar wielokąta, statystyki strefowe ...), które mogą być przydatne do dalszego przetwarzania.

Alternatywy 2b) W zależności od potrzeb możesz wybrać lokalizacje otaczających wielokątów w odległości 200 m, a następnie zidentyfikować niektóre rodzaje „obudów”, ale nie tak ściśle jak w 2).

Biorąc pod uwagę „przypadek labiryntu”, może to pomóc: ocenić, jak długo trzeba „uciekać” z lokalizacji.

• Możesz już wykluczyć z analizy punkty, które są w pełni uwzględnione lub całkowicie bezpłatne

• następnie przekształcasz swoje przeszkody w raster i ustawiasz wartości na NoData, gdzie masz wielokąt, i na wielkość komórki w metrach, gdzie nie masz (spowoduje to raster kosztów).

• po trzecie, możesz obliczyć odległość kosztów za pomocą nowo wygenerowanego rastra kosztów

• na koniec używasz statystyki strefowej jako tabeli opartej na granicach bufora przekonwertowanych na raster (tworząc pierścień). Jeśli możesz uciec we wszystkich kierunkach, minimum powinno wynosić około 200 (w zależności od wielkości komórki analizy). Ale jeśli jesteś w labiryncie, maksimum będzie większe niż 200. Zatem maksimum statystyk strefowych minus 200 będzie ciągłą wartością wskazującą, jak trudno jest „uciec”.