Wie erstelle ich einen orientierten Puffer mit arcpy?

Ich möchte mit arcpy für jedes Polygon in meinem Shapefile einen orientierten Puffer erstellen. Mit orientiert meine ich, dass ich zwei Winkel a1 und a2 habe, die die Richtung des Puffers einschränken. Dies ist in der folgenden Grafik dargestellt:

Irgendwelche Ideen?

Zusammenfassung

Diese Antwort stellt die Frage in einen größeren Kontext, beschreibt einen effizienten Algorithmus, der auf die Shapefile-Darstellung von Features (als "Vektoren" oder "Linestrings" von Punkten) anwendbar ist, zeigt einige Beispiele für ihre Anwendung und gibt Arbeitscode für die Verwendung oder Portierung an eine GIS-Umgebung.

Hintergrund

Dies ist ein Beispiel für eine morphologische Erweiterung. Im Allgemeinen "verbreitet" eine Erweiterung die Punkte einer Region in ihre Nachbarschaften; Die Sammlung von Punkten, an denen sie landen, ist die "Erweiterung". Die Anwendungen in GIS sind zahlreich: Modellierung der Ausbreitung von Feuer, der Bewegung von Zivilisationen, der Ausbreitung von Pflanzen und vielem mehr.

Mathematisch und in sehr großer (aber nützlicher) Allgemeinheit breitet eine Dilatation eine Reihe von Punkten in einer Riemannschen Mannigfaltigkeit aus (z. B. eine Ebene, eine Kugel oder ein Ellipsoid). Die Ausbreitung wird durch eine Teilmenge des Tangentenbündels an diesen Punkten festgelegt. Dies bedeutet, dass an jedem der Punkte eine Reihe von Vektoren (Richtungen und Entfernungen) angegeben ist (ich nenne dies eine "Nachbarschaft"); Jeder dieser Vektoren beschreibt einen geodätischen Pfad, der an seinem Basispunkt beginnt. Der Basispunkt wird bis zu den Enden all dieser Pfade "verteilt". (Für die viel eingeschränktere Definition von "Dilatation", die üblicherweise in der Bildverarbeitung verwendet wird, siehe den Wikipedia-Artikel . Die Ausbreitungsfunktion wird als Exponentialkarte bezeichnet in Differentialgeometrie.)

Das "Puffern" eines Merkmals ist eines der einfachsten Beispiele für eine solche Erweiterung: Um jeden Punkt des Merkmals wird (zumindest konzeptionell) eine Scheibe mit konstantem Radius (der Pufferradius) erstellt. Die Vereinigung dieser Festplatten ist der Puffer.

Diese Frage erfordert die Berechnung einer etwas komplizierteren Dilatation, bei der die Ausbreitung nur innerhalb eines bestimmten Winkelbereichs (dh innerhalb eines kreisförmigen Sektors) erfolgen darf. Dies ist nur für Merkmale sinnvoll, die keine merklich gekrümmte Oberfläche einschließen (z. B. kleine Merkmale auf der Kugel oder dem Ellipsoid oder Merkmale in der Ebene). Wenn wir in der Ebene arbeiten, ist es auch sinnvoll, alle Sektoren in die gleiche Richtung auszurichten. (Wenn wir jedoch die Ausbreitung von Feuer durch Wind modellieren würden, möchten wir, dass die Sektoren auf den Wind ausgerichtet sind und ihre Größe auch mit der Windgeschwindigkeit variieren kann: Dies ist eine Motivation für die allgemeine Definition der Dilatation, die ich gegeben habe. ) (Auf gekrümmten Flächen wie einem Ellipsoid ist es im Allgemeinen unmöglich, alle Sektoren in der "gleichen" Richtung auszurichten.)

Unter folgenden Umständen ist die Dilatation relativ einfach zu berechnen:

• Das Feature befindet sich in der Ebene (das heißt, wir erweitern eine Karte des Features und hoffentlich ist die Karte ziemlich genau).

• Die Dilatation ist konstant : Die Ausbreitung an jedem Punkt des Merkmals erfolgt in kongruenten Nachbarschaften mit derselben Ausrichtung.

• Diese gemeinsame Nachbarschaft ist konvex. Die Konvexität vereinfacht und beschleunigt die Berechnungen erheblich.

Diese Frage passt zu solchen speziellen Umständen: Sie fordert die Erweiterung beliebiger Polygone durch kreisförmige Sektoren, deren Ursprung (die Zentren der Scheiben, von denen sie stammen) an den Basispunkten liegt. Vorausgesetzt, diese Sektoren erstrecken sich nicht über mehr als 180 Grad, sind sie konvex. (Größere Sektoren können immer in zwei konvexe Sektoren aufgeteilt werden. Die Vereinigung der beiden kleineren Dilatationen ergibt das gewünschte Ergebnis.)

Implementierung

Weil wir die euklidische Berechnungen durchführen - dabei in der Ebene der Ausbreitung - können wir einen Punkt aufzuweiten lediglich durch die Übersetzung der Erweiterung Nachbarschaft zu diesem Punkt. (Um dies tun zu können, braucht die Nachbarschaft einen Ursprungdas entspricht dem Basispunkt. Zum Beispiel ist der Ursprung der Sektoren in dieser Frage der Mittelpunkt des Kreises, aus dem sie gebildet werden. Dieser Ursprung liegt zufällig an der Grenze des Sektors. Bei der Standard-GIS-Pufferoperation ist die Nachbarschaft ein Kreis mit dem Ursprung in der Mitte. jetzt liegt der Ursprung im Inneren des Kreises. Die Auswahl eines Ursprungs ist rechnerisch keine große Sache, da eine Änderung des Ursprungs lediglich die gesamte Dilatation verschiebt, aber es kann eine große Sache in Bezug auf die Modellierung natürlicher Phänomene sein. Die sectorFunktion im folgenden Code zeigt, wie ein Ursprung angegeben werden kann.)

Das Erweitern eines Liniensegments kann schwierig sein, aber für eine konvexe Nachbarschaft können wir die Erweiterung als Vereinigung der Erweiterungen der beiden Endpunkte zusammen mit einem sorgfältig ausgewählten Parallelogramm erstellen. (Im Interesse des Weltraums werde ich nicht innehalten, um solche mathematischen Behauptungen zu beweisen, sondern die Leser ermutigen, ihre eigenen Beweise zu versuchen, da dies eine aufschlussreiche Übung ist.) Hier ist eine Illustration mit drei Sektoren (in Pink dargestellt). Sie haben Einheitsradien und ihre Winkel sind in den Titeln angegeben. Das Liniensegment selbst hat die Länge 2, ist horizontal und wird schwarz angezeigt:

Die Parallelogramme werden gefunden, indem die rosa Punkte, die so weit wie möglich vom Segment entfernt sind, nur in vertikaler Richtung lokalisiert werden . Dies ergibt zwei untere Punkte und zwei obere Punkte entlang Linien, die parallel zum Segment sind. Wir müssen nur die vier Punkte zu einem Parallelogramm zusammenfügen (blau dargestellt). Beachten Sie rechts, wie sinnvoll dies ist, auch wenn der Sektor selbst nur ein Liniensegment (und kein echtes Polygon) ist: Dort wurde jeder Punkt auf dem Segment für eine Entfernung von 171 Grad östlich von Norden verschoben von 0 bis 1. Die Menge dieser Endpunkte ist das gezeigte Parallelogramm. Die Details dieser Berechnung werden in der im folgenden Code bufferdefinierten Funktion angezeigt dilate.edges.

Um eine Polylinie zu erweitern , bilden wir die Vereinigung der Erweiterungen der Punkte und Segmente, die sie bilden. Die letzten beiden Zeilen dilate.edgesführen diese Schleife aus.

Um ein Polygon zu erweitern, muss das Innere des Polygons zusammen mit der Erweiterung seiner Grenze berücksichtigt werden . (Diese Behauptung macht einige Annahmen über die Dilatationsnachbarschaft. Eine ist, dass alle Nachbarschaften den Punkt (0,0) enthalten, der garantiert, dass das Polygon in seiner Dilatation enthalten ist. Bei variablen Nachbarschaften wird auch die Dilatation eines Innenraums angenommen Der Punkt des Polygons erstreckt sich nicht über die Erweiterung der Grenzpunkte hinaus. Dies ist bei konstanten Nachbarschaften der Fall.)

Schauen wir uns einige Beispiele an, wie dies funktioniert, zuerst mit einem Nonagon (ausgewählt, um Details zu enthüllen) und dann mit einem Kreis (ausgewählt, um der Abbildung in der Frage zu entsprechen). Die Beispiele werden weiterhin dieselben drei Nachbarschaften verwenden, jedoch auf einen Radius von 1/3 geschrumpft.

In dieser Abbildung ist das Innere des Polygons grau, Punktdilatationen (Sektoren) sind rosa und Randdilatationen (Parallelogramme) sind blau.

Der "Kreis" ist eigentlich nur ein 60-Gon, aber er nähert sich gut einem Kreis.

Performance

Wenn das Basismerkmal durch N Punkte und die Dilatationsnachbarschaft durch M Punkte dargestellt wird, erfordert dieser Algorithmus O (N M) Aufwand. Daraufhin muss das Durcheinander von Scheitelpunkten und Kanten in der Vereinigung vereinfacht werden, was O (N M log (N M)) Aufwand erfordern kann: Das muss das GIS tun. das sollten wir nicht programmieren müssen.

Der Rechenaufwand für konvexe Basismerkmale könnte auf O (M + N) verbessert werden (da Sie herausfinden können, wie Sie die neue Grenze umgehen, indem Sie die Scheitelpunktlisten, die die Grenzen der beiden ursprünglichen Formen beschreiben, entsprechend zusammenführen). Dies würde auch keine anschließende Reinigung erfordern.

Wenn sich die Größe und / oder Ausrichtung der Dilatationsumgebung langsam ändert, während Sie sich um das Basismerkmal bewegen, kann die Dilatation der Kante von der konvexen Hülle der Vereinigung der Dilatationen ihrer Endpunkte aus genau angenähert werden. Wenn die beiden Dilatationsnachbarschaften M1- und M2-Punkte haben, kann dies mit O (M1 + M2) -Anstrengung unter Verwendung eines in Shamos & Preparata, Computational Geometry beschriebenen Algorithmus gefunden werden . Wenn also K = M1 + M2 + ... + M (N) die Gesamtzahl der Eckpunkte in den N Dilatationsvierteln ist, können wir die Dilatation in O (K * log (K)) Zeit berechnen.

Warum sollten wir eine solche Verallgemeinerung in Angriff nehmen wollen, wenn wir nur einen einfachen Puffer wollen? Bei großen Merkmalen auf der Erde kann eine Dilatationsumgebung (z. B. eine Scheibe) mit einer in der Realität konstanten Größe auf der Karte, auf der diese Berechnungen durchgeführt werden, eine unterschiedliche Größe haben. Auf diese Weise erhalten wir eine Möglichkeit, Berechnungen durchzuführen, die für das Ellipsoid genau sind, während wir weiterhin alle Vorteile der euklidischen Geometrie nutzen können.

Code

Die Beispiele wurden mit diesem RPrototyp erstellt, der problemlos in Ihre Lieblingssprache (Python, C ++ usw.) portiert werden kann. In der Struktur entspricht es der in dieser Antwort angegebenen Analyse und bedarf daher keiner gesonderten Erläuterung. Kommentare verdeutlichen einige Details.

(Es kann zu Notiz interessant sein , dass trigonometrische Berechnungen verwendet werden , nur die Beispielfunktionen zu schaffen - die regelmäßige Polygone sind -. Und die Sektoren Kein Teil der Dilatation Berechnungen erfordert jede Trigonometrie.)

#
# Dilate the vertices of a polygon/polyline by a shape.
#
dilate.points <- function(p, q) {
  # Translate a copy of `q` to each vertex of `p`, resulting in a list of polygons.
  pieces <- apply(p, 1, function(x) list(t(t(q)+x)))
  lapply(pieces, function(z) z[[1]]) # Convert to a list of matrices
}
#
# Dilate the edges of a polygon/polyline `p` by a shape `q`. 
# `p` must have at least two rows.
#
dilate.edges <- function(p, q) {
  i <- matrix(c(0,-1,1,0), 2, 2)       # 90 degree rotation
  e <- apply(rbind(p, p[1,]), 2, diff) # Direction vectors of the edges
  # Dilate a single edge from `x` to `x+v` into a parallelogram
  # bounded by parts of the dilation shape that are at extreme distances
  # from the edge.
  buffer <- function(x, v) {
    y <- q %*% i %*% v # Signed distances orthogonal to the edge
    k <- which.min(y)  # Find smallest distance, then the largest *after* it
    l <- (which.max(c(y[-(1:k)], y[1:k])) + k-1) %% length(y)[1] + 1
    list(rbind(x+q[k,], x+v+q[k,], x+v+q[l,], x+q[l,])) # A parallelogram
  }
  # Apply `buffer` to every edge.
  quads <- apply(cbind(p, e), 1, function(x) buffer(x[1:2], x[3:4]))
  lapply(quads, function(z) z[[1]]) # Convert to a list of matrices
}
#----------------------- (This ends the dilation code.) --------------------------#
#
# Display a polygon and its point and edge dilations.
# NB: In practice we would submit the polygon, its point dilations, and edge 
#     dilations to the GIS to create and simplify their union, producing a single
#     polygon.  We keep the three parts separate here in order to illustrate how
#     that polygon is constructed.
#
display <- function(p, d.points, d.edges, ...) {
  # Create a plotting region covering the extent of the dilated figure.
  x <- c(p[,1], unlist(lapply(c(d.points, d.edges), function(x) x[,1])))
  y <- c(p[,2], unlist(lapply(c(d.points, d.edges), function(x) x[,2])))
  plot(c(min(x),max(x)), c(min(y),max(y)), type="n", asp=1, xlab="x", ylab="y", ...)
  # The polygon itself.
  polygon(p, density=-1, col="#00000040")
  # The dilated points and edges.
  plot.list <- function(l, c) lapply(l, function(p) 
                  polygon(p, density=-1, col=c, border="#00000040"))
  plot.list(d.points, "#ff000020")
  plot.list(d.edges, "#0000ff20")
  invisible(NULL) # Doesn't return anything
}
#
# Create a sector of a circle.
# `n` is the number of vertices to use for approximating its outer arc.
#
sector <- function(radius, arg1, arg2, n=1, origin=c(0,0)) {
  t(cbind(origin, radius*sapply(seq(arg1, arg2, length.out=n), 
                  function(a) c(cos(a), sin(a)))))
}
#
# Create a polygon represented as an array of rows.
#
n.vertices <- 60 # Inscribes an `n.vertices`-gon in the unit circle.
angles <- seq(2*pi, 0, length.out=n.vertices+1)
angles <- angles[-(n.vertices+1)]
polygon.the <- cbind(cos(angles), sin(angles))
if (n.vertices==1) polygon.the <- rbind(polygon.the, polygon.the)
#
# Dilate the polygon in various ways to illustrate.
#
system.time({
  radius <- 1/3
  par(mfrow=c(1,3))
  q <- sector(radius, pi/12, 2*pi/3, n=120)
  d.points <- dilate.points(polygon.the, q)
  d.edges <- dilate.edges(polygon.the, q)
  display(polygon.the, d.points, d.edges, main="-30 to 75 degrees")

  q <- sector(radius, pi/3, 4*pi/3, n=180)
  d.points <- dilate.points(polygon.the, q)
  d.edges <- dilate.edges(polygon.the, q)
  display(polygon.the, d.points, d.edges, main="-150 to 30 degrees")

  q <- sector(radius, -9/20*pi, -9/20*pi)
  d.points <- dilate.points(polygon.the, q)
  d.edges <- dilate.edges(polygon.the, q)
  display(polygon.the, d.points, d.edges, main="171 degrees")
})

Die Rechenzeit für dieses Beispiel (aus der letzten Abbildung) mit N = 60 und M = 121 (links), M = 181 (Mitte) und M = 2 (rechts) betrug eine Viertelsekunde. Das meiste davon war jedoch für die Anzeige. Normalerweise verarbeitet dieser RCode ungefähr N M = 1,5 Millionen pro Sekunde (es dauert nur etwa 0,002 Sekunden, um alle gezeigten Beispielberechnungen durchzuführen). Das Erscheinungsbild des Produkts M N impliziert jedoch, dass die Erweiterung vieler Figuren oder komplizierter Figuren über eine detaillierte Nachbarschaft beträchtliche Zeit in Anspruch nehmen kann. Seien Sie also vorsichtig! Vergleichen Sie das Timing für kleinere Probleme, bevor Sie ein großes Problem angehen. Unter solchen Umständen könnte man nach einer rasterbasierten Lösung suchen (die viel einfacher zu implementieren ist und im Wesentlichen nur eine einzige Nachbarschaftsberechnung erfordert).

Das ist ziemlich weit gefasst, aber Sie könnten:

1. Puffern Sie das ursprüngliche Polygon
2. Finden Sie den Ursprungspunkt der "orientierten" Strahlen, die an der Polygongrenze erzeugt werden sollen (eine Art Tangentialpunkt?)
3. Erstellen / Erweitern einer Linie von diesem Punkt bis zu einer Entfernung jenseits des Puffers unter Verwendung des in den Kommentaren zur Frage beschriebenen Winkels.
4. Schneiden Sie diese Linie mit dem Puffer und dem ursprünglichen Polygon. Dies könnte wahrscheinlich gleichzeitig mit 3) mit den richtigen Argumenten für Extend geschehen.
5. Extrahieren Sie das neue Polygon "Orientierter Puffer" aus der resultierenden Menge von Polygonen