Feststellen, ob der Punkt mithilfe der Rasterverarbeitung umgeben ist

Ich versuche, einen derzeit äußerst umständlichen Vektor / Python-Prozess für ein Naturgefahrenmodell zu verbessern. Im Moment haben wir ein langes Skript, das Entfernungs- / Peilungslinien von einem bestimmten Punkt generiert, um Folgendes zu bestimmen:

1. die Art des Polygons, das es schneidet (z. B. Wald, Gras, Sumpf usw.)
2. der Abstand zu diesem Polygon
3. Wie viele dieser Linien schneiden Polygone, um festzustellen, wie "umgeben" sie sind.

Es geht um viel mehr, aber das ist das Wesentliche. Ich versuche einen Weg zu finden, um dies zu verbessern, und bin derzeit bei Teil 3 ratlos. Die Idee ist, festzustellen, ob ein Punkt innerhalb von beispielsweise 200 m vollständig von Polygonen umgeben ist

In meinem angehängten Bild möchte ich, dass Punkt A einem höheren Risiko ausgesetzt ist als Punkt B, da er vollständig von meinen Polygonen umgeben ist. Dies wird für ungefähr 13 Millionen Punkte wiederholt, ist also keine kleine Aufgabe, und ich hätte lieber eine Oberfläche, aus der ich Werte ableiten kann, als unser Skript auszuführen. Ich denke, es muss eine Variation von Hydrologie-Tools oder Kostenpfaden geben, um dies zu tun, aber ich kann mich anscheinend nicht darum kümmern.

Wie könnte ich das machen?

Es gibt eine Kostenpfadlösung, die Sie jedoch selbst codieren müssen. So könnte es aussehen, wenn es auf jeden Punkt im Bild in der Frage angewendet wird (etwas vergröbert, um die Berechnungen zu beschleunigen):

Die schwarzen Zellen sind Teile der umgebenden Polygone. Die Farben, die von hellorange (kurz) bis blau (lang) reichen, zeigen die maximale Entfernung (bis zu maximal 50 Zellen) an, die durch Sichtlinienüberquerung erreicht werden kann, ohne die Polygonzellen abzufangen. (Jede Zelle außerhalb der Ausdehnung dieses Bildes wird als Teil der Polygone behandelt.)

Lassen Sie uns einen effizienten Weg diskutieren, dies mithilfe einer Rasterdarstellung der Daten zu tun. In dieser Darstellung haben alle "umgebenden" polygonalen Zellen beispielsweise Werte ungleich Null, und jede Zelle, die "durchschaut" werden kann, hat einen Wert von Null.

Schritt 1: Vorberechnung einer Nachbarschaftsdatenstruktur

Sie müssen zuerst entscheiden, was es für eine Zelle bedeutet, eine andere zu blockieren. Eine der fairsten Regeln, die ich finden kann, lautet: Verwenden Sie Integralkoordinaten für Zeilen und Spalten (und nehmen Sie quadratische Zellen an), und überlegen Sie, welche Zellen die Zelle (i, j) aus der Ansicht am Ursprung (0,0) blockieren könnten. Ich nominiere die Zelle (i ', j'), die dem Liniensegment am nächsten liegt, das (i, j) mit (0,0) verbindet, und zwar zwischen allen Zellen, deren Koordinaten sich von i und j um höchstens 1 unterscheiden. Weil dies nicht immer der Fall ist ergeben eine eindeutige Lösung (zum Beispiel mit (i, j) = (1,2) funktionieren sowohl (0,1) als auch (1,1) gleich gut), einige Mittel zum Auflösen von Bindungen sind erforderlich. Es wäre schön, wenn diese Auflösung von Bindungen die Symmetrien kreisförmiger Nachbarschaften in Gittern berücksichtigen würde: Wenn Sie entweder die Koordinaten negieren oder die Koordinaten wechseln, bleiben diese Nachbarschaften erhalten. Daher können wir entscheiden, welche Zellen blockieren (i,

Zur Veranschaulichung dieser Regel wird der folgende Prototypcode geschrieben R. Dieser Code gibt eine Datenstruktur zurück, die zur Bestimmung der "Umgebung" beliebiger Zellen in einem Raster geeignet ist.

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

Der Wert von screen(12)wurde verwendet, um diese Darstellung dieser Screening-Beziehung zu erstellen: Pfeile zeigen von den Zellen zu denen, die sie sofort screenen. Die Farbtöne sind proportional zur Entfernung zum Ursprung, der sich in der Mitte dieses Viertels befindet:

Diese Berechnung ist schnell und muss für eine bestimmte Nachbarschaft nur einmal durchgeführt werden. Wenn Sie beispielsweise 200 m in einem Raster mit 5 m Zellen betrachten, beträgt die Nachbarschaftsgröße 200/5 = 40 Einheiten.

Schritt 2: Anwenden der Berechnung auf ausgewählte Punkte

Der Rest ist unkompliziert: Um festzustellen, ob eine Zelle an (x, y) (in Zeilen- und Spaltenkoordinaten) in Bezug auf diese Nachbarschaftsdatenstruktur "umgeben" ist, führen Sie den Test rekursiv aus, beginnend mit einem Versatz von (i, j). = (0,0) (der Nachbarschaftsursprung). Wenn der Wert im Polygonraster bei (x, y) + (i, j) ungleich Null ist, ist die Sichtbarkeit dort blockiert. Andernfalls müssen wir alle Offsets berücksichtigen, die beim Offset (i, j) blockiert worden sein könnten (die in O (1) -Zeit unter Verwendung der von zurückgegebenen Datenstruktur gefunden werden screen). Wenn keine blockiert sind, haben wir den Umfang erreicht und schließen daraus, dass (x, y) nicht umgeben ist. Daher stoppen wir die Berechnung (und machen uns nicht die Mühe, verbleibende Punkte in der Nachbarschaft zu untersuchen).

Wir können noch nützlichere Informationen sammeln, indem wir die weiteste Sichtlinienentfernung verfolgen, die während des Algorithmus erreicht wurde. Wenn dies kleiner als der gewünschte Radius ist, ist die Zelle umgeben; sonst ist es nicht.

Hier ist ein RPrototyp dieses Algorithmus. Es ist länger als es scheint, da Res die (einfache) Stapelstruktur, die zum Implementieren der Rekursion erforderlich ist, nicht nativ unterstützt, sodass auch ein Stapel codiert werden muss. Der eigentliche Algorithmus beginnt ungefähr zwei Drittel des Weges und benötigt nur etwa ein Dutzend Zeilen. (Und die Hälfte von ihnen kümmert sich lediglich um die Situation am Rand des Rasters und sucht nach Indizes außerhalb des Bereichs in der Nachbarschaft. Dies könnte effizienter gestaltet werden, indem das Polygonraster einfach um kZeilen und Spalten um seinen Umfang erweitert wird, wodurch alle eliminiert werden Notwendigkeit einer Überprüfung des Indexbereichs auf Kosten von etwas mehr RAM, um das Polygongitter zu halten.)

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

In diesem Beispiel sind polygonale Zellen schwarz. Farben geben den maximalen Sichtlinienabstand (bis zu 50 Zellen) für nicht polygonale Zellen an, der von hellorange für kurze Entfernungen bis dunkelblau für die längsten Entfernungen reicht. (Die Zellen sind eine Einheit breit und hoch.) Die sichtbar sichtbaren Streifen werden durch die kleinen Polygon- "Inseln" in der Mitte des "Flusses" erzeugt: Jede blockiert eine lange Reihe anderer Zellen.

Analyse des Algorithmus

Die Stapelstruktur implementiert eine Tiefensuche des Nachbarschaftssichtbarkeitsgraphen, um zu beweisen, dass eine Zelle nicht umgeben ist. Wenn Zellen weit von einem Polygon entfernt sind, müssen bei dieser Suche nur O (k) -Zellen auf eine kreisförmige Nachbarschaft mit Radius-k untersucht werden. Die schlimmsten Fälle treten auf, wenn es eine kleine Anzahl verstreuter Polygonzellen in der Nachbarschaft gibt, aber dennoch die Grenze der Nachbarschaft nicht ganz erreicht werden kann: Diese erfordern die Inspektion fast aller Zellen in jeder Nachbarschaft, was ein O (k ^ 2) ist. Betrieb.

Das folgende Verhalten ist typisch für das, was angetroffen wird. Bei kleinen Werten von k sind die meisten nicht polygonalen Zellen offensichtlich nicht umgeben, sofern die Polygone den größten Teil des Gitters nicht ausfüllen, und der Algorithmus skaliert wie O (k). Bei Zwischenwerten sieht die Skalierung wie O (k ^ 2) aus. Wenn k wirklich groß wird, werden die meisten Zellen umgeben sein und diese Tatsache kann lange vor der Inspektion der gesamten Nachbarschaft bestimmt werden: Der Rechenaufwand des Algorithmus erreicht dadurch eine praktische Grenze. Diese Grenze wird erreicht, wenn sich der Nachbarschaftsradius dem Durchmesser der größten verbundenen nicht-polygonalen Bereiche im Gitter nähert.

Als Beispiel verwende ich die countingim Prototyp von codierte Option screen, um die Anzahl der in jedem Aufruf verwendeten Stapeloperationen zurückzugeben. Dies misst den Rechenaufwand. Das folgende Diagramm zeigt die mittlere Anzahl von Stack-Ops als Funktion des Nachbarschaftsradius. Es zeigt das vorhergesagte Verhalten.

Wir können dies verwenden, um die Berechnung zu schätzen, die zur Bewertung von 13 Millionen Punkten in einem Raster erforderlich ist. Angenommen, eine Nachbarschaft von k = 200/5 = 40 wird verwendet. Dann werden durchschnittlich einige hundert Stapeloperationen benötigt (abhängig von der Komplexität des Polygongitters und der Position der 13 Millionen Punkte relativ zu den Polygonen), was bedeutet, dass in einer effizienten kompilierten Sprache höchstens einige tausend einfache numerische Operationen ausgeführt werden wird benötigt (addieren, multiplizieren, lesen, schreiben, versetzen usw.). Die meisten PCs werden in der Lage sein, die Umgebung von etwa einer Million Punkten bei dieser Rate zu bewerten. (DasRDie Implementierung ist viel, viel langsamer als diese, da sie bei dieser Art von Algorithmus schlecht ist, weshalb sie nur als Prototyp betrachtet werden kann.) Dementsprechend können wir hoffen, dass eine effiziente Implementierung in einer einigermaßen effizienten und geeigneten Sprache - C ++ - möglich ist und Python kommen in den Sinn - könnte die Auswertung von 13 Millionen Punkten in einer Minute oder weniger abschließen, vorausgesetzt, das gesamte Polygongitter befindet sich im RAM.

Wenn ein Raster zu groß ist, um in den Arbeitsspeicher zu passen, kann dieses Verfahren auf gekachelte Teile des Rasters angewendet werden. Sie müssen sich nur durch kZeilen und Spalten überlappen . Nehmen Sie die Maxima an den Überlappungen, wenn Sie die Ergebnisse mosaikieren.

Andere Anwendungen

Das "Holen" eines Gewässers hängt eng mit der "Umgebung" seiner Punkte zusammen. Wenn wir einen Nachbarschaftsradius verwenden, der gleich oder größer als der Durchmesser des Wasserkörpers ist, erstellen wir an jedem Punkt im Wasserkörper ein Gitter des (nicht gerichteten) Abrufs. Durch Verwendung eines kleineren Nachbarschaftsradius erhalten wir zumindest eine Untergrenze für den Abruf an allen Punkten mit dem höchsten Abruf, was in einigen Anwendungen gut genug sein kann (und den Rechenaufwand erheblich reduzieren kann). Eine Variante dieses Algorithmus, die die "gescreent von" -Relation auf bestimmte Richtungen begrenzt, wäre eine Möglichkeit, den Abruf in diesen Richtungen effizient zu berechnen. Beachten Sie, dass für solche Varianten der Code geändert werden muss screen. Der Code für panvisibilityändert sich überhaupt nicht.

Ich kann definitiv sehen, wie man dies mit einer Rasterlösung machen möchte, aber selbst bei einer reduzierten Anzahl von Punkten würde ich ein sehr großes / hochauflösendes und daher schwer zu verarbeitendes Gitter oder eine Reihe von Gittern erwarten. Angesichts dessen frage ich mich, ob das Ausnutzen der Topologie in einer GDB effizienter sein könnte. Sie könnten alle internen Hohlräume mit so etwas wie finden:

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')`

Sie können dann topoErrorsVoidPolysin Ihrem normalen Muster von Intersect_analysis()oder was auch immer arbeiten. Möglicherweise müssen Sie mit dem Extrahieren der Polys von Interesse herumspielen topoErrorsVoidPolys. @whuber hat eine Reihe ziemlich ausgezeichneter Beiträge zu solchen Dingen an anderer Stelle hier auf gis.stackexchange.com.

Wenn du wirklich Raster machen willst ... würde ich etwas in der Art dieses Pseudocodes tun (nicht zusammenzucken, nur weil es offensichtlich ist, dass ich ein AML-Rückfall bin !: P)

1. rastern Sie Punkte ("pts_g") und polys ("polys_g" (
2. Hohlräume = Regionsgruppe (con (isnull (polys_g), 1))
3. Möglicherweise müssen Sie etwas tun, um Hohlräume zu verfeinern und unerwünschte externe Polygone / offene Universumsbereiche zu beseitigen
4. pts_surrounded = con (Hohlräume, pts_g)

Machen Sie das einfach nach, also müssen Sie es vielleicht verfeinern.

Wenn Sie einen Schwellenabstandswert verwenden (hier sprechen Sie von 200 m), ist die beste Lösung die Vektoranalyse:

1) Erstellen Sie einen 200-m-Puffer um jeden Punkt (in der Abbildung schwarz).

2) Verwenden Sie das Schnittwerkzeug (Analyse) zwischen dem Puffer und den Polygonen (in der Abbildung blau). Es sieht besser aus, wenn Sie dies zwischen den Grenzen der umgebenden Polygone und dem Puffer tun, aber das Endergebnis ist das gleiche.

3) Verwenden Sie die Funktion zum Polygonieren (Verwaltung), um Polygone zu erstellen, bei denen Ihre Punkte vollständig umgeben sind (in der Abbildung rot).

4) Wählen Sie Ebenen nach Standort (Verwaltung) oder räumlicher Verknüpfung (Analyse) aus, um die umgebenden Punkte zu identifizieren. Durch die Verwendung der räumlichen Verknüpfung erhalten Sie Informationen zum Einbettungspolygon (Bereich des Polygons, Zonenstatistik ...), die für die weitere Verarbeitung hilfreich sein können.

Alternativen 2b) Abhängig von Ihren Anforderungen können Sie die umgebenden Polygone in einem Abstand von 200 m nach Standort auswählen. Anschließend können Sie einige Arten von "Gehäusen" identifizieren, die jedoch nicht so streng sind wie in 2).

In Anbetracht des "Labyrinthfalls" könnte dies helfen: Bewerten Sie, wie lange es dauert, um vom Ort "zu entkommen".

• Sie können bereits Punkte, die vollständig enthalten oder vollständig frei sind, von der Analyse ausschließen

• Dann konvertieren Sie Ihre Hindernisse in ein Raster und setzen die Werte auf NoData, wo Sie ein Polygon haben, und auf die Zellengröße in Metern, wo Sie dies nicht tun (dies macht Ihr Kostenraster).

• Drittens können Sie die Kostenentfernung mit dem neu generierten Kostenraster berechnen

• Schließlich verwenden Sie eine Zonenstatistik als Tabelle basierend auf den in Raster konvertierten Puffergrenzen (Bildung eines Rings). Wenn Sie in alle Richtungen entkommen können, sollte das Minimum ungefähr 200 betragen (abhängig von der Zellengröße Ihrer Analyse). Wenn Sie sich jedoch in einem Labyrinth befinden, ist das Maximum größer als 200. Das Maximum der Zonenstatistik minus 200 ist also ein kontinuierlicher Wert, der angibt, wie schwierig es ist, "zu entkommen".