Berechnung der maximalen Entfernung innerhalb eines Polygons in x-Richtung (Ost-West-Richtung) in PostGIS?

Ich interessiere mich für die maximale Breite eines Polygons, zB eines Sees in Ost-West-Richtung. Begrenzungsrahmen helfen nur bei einfachen Polygonen, nicht aber bei komplexen konkaven Polygonen.

Dies erfordert wahrscheinlich einige Skripte in jeder GIS-Plattform.

Die effizienteste Methode (asymptotisch) ist eine vertikale Linienabtastung: Sie erfordert das Sortieren der Kanten nach ihren minimalen y-Koordinaten und das anschließende Verarbeiten der Kanten von unten (minimales y) nach oben (maximales y) für ein O (e * log ( e)) Algorithmus, wenn e Kanten betroffen sind.

Obwohl das Verfahren einfach ist, ist es überraschend schwierig, es in allen Fällen richtig zu machen. Polygone können böse sein: Sie können baumeln, splittern, Löcher aufweisen, voneinander getrennt sein, doppelte Scheitelpunkte aufweisen, Scheitelpunkte entlang gerader Linien verlaufen und ungelöste Grenzen zwischen zwei benachbarten Komponenten aufweisen. Hier ist ein Beispiel mit vielen dieser Eigenschaften (und mehr):

Wir werden speziell das horizontale Segment maximaler Länge suchen, das vollständig innerhalb des Verschlusses des Polygons liegt. Dies beseitigt beispielsweise das Baumeln zwischen x = 20 und x = 40, das von dem Loch zwischen x = 10 und x = 25 ausgeht. Es ist dann einfach zu zeigen, dass mindestens eines der horizontalen Segmente maximaler Länge mindestens einen Scheitelpunkt schneidet. (Wenn es Lösungen gibt keine Ecken schneiden, werden sie in das Innere einiger Parallelogramm an der Ober- und Unterseite von Lösungen begrenzt liegen , die nicht schneiden mindestens eine Ecke. Das gibt uns ein Mittel zu finden , alle Lösungen.)

Dementsprechend muss der Linien-Sweep mit den niedrigsten Scheitelpunkten beginnen und sich dann nach oben bewegen (d. H. Zu höheren y-Werten), um an jedem Scheitelpunkt anzuhalten. An jeder Haltestelle finden wir neue Kanten, die von dieser Höhe nach oben ragen. Eliminieren Sie alle Kanten, die in dieser Höhe von unten enden (dies ist eine der Schlüsselideen: Es vereinfacht den Algorithmus und eliminiert die Hälfte der potenziellen Verarbeitung). und sorgfältig alle Kanten bearbeiten, die vollständig in einer konstanten Höhe liegen (die horizontalen Kanten).

Betrachten Sie beispielsweise den Zustand, in dem ein Niveau von y = 10 erreicht ist. Von links nach rechts finden wir folgende Kanten:

      x.min x.max y.min y.max
 [1,]    10     0     0    30
 [2,]    10    24    10    20
 [3,]    20    24    10    20
 [4,]    20    40    10    10
 [5,]    40    20    10    10
 [6,]    60     0     5    30
 [7,]    60    60     5    30
 [8,]    60    70     5    20
 [9,]    60    70     5    15
[10,]    90   100    10    40

In dieser Tabelle sind (x.min, y.min) Koordinaten des unteren Endpunkts der Kante und (x.max, y.max) Koordinaten ihres oberen Endpunkts. Auf dieser Ebene (y = 10) wird die erste Kante in ihrem Inneren abgefangen, die zweite Kante unten und so weiter. Einige Kanten, die auf dieser Ebene enden, wie z. B. (10,0) bis (10,10), sind nicht in der Liste enthalten.

Um festzustellen, wo sich die inneren und äußeren Punkte befinden, stellen Sie sich vor, Sie beginnen ganz links - natürlich außerhalb des Polygons - und bewegen sich horizontal nach rechts. Jedes Mal, wenn wir eine Kante überqueren, die nicht horizontal ist , wechseln wir abwechselnd von außen nach innen und zurück. (Dies ist eine weitere Schlüsselidee.) Es wird jedoch bestimmt, dass alle Punkte innerhalb einer horizontalen Kante innerhalb des Polygons liegen, egal was passiert. (Der Abschluss eines Polygons schließt immer seine Kanten ein.)

Im folgenden Beispiel sehen Sie die sortierte Liste der x-Koordinaten, bei denen nicht horizontale Kanten an der y = 10-Linie beginnen oder diese kreuzen:

x.array    6.7 10 20 48 60 63.3 65 90
interior     1  0  1  0  1    0  1  0

(Beachten Sie, dass x = 40 nicht in dieser Liste enthalten ist.) Die Werte des interiorArrays markieren die linken Endpunkte der inneren Segmente: 1 bezeichnet ein inneres Intervall, 0 ein äußeres Intervall. Somit gibt die erste 1 an, dass das Intervall von x = 6,7 bis x = 10 innerhalb des Polygons liegt. Die nächste 0 gibt an, dass das Intervall von x = 10 bis x = 20 außerhalb des Polygons liegt. Und so geht es weiter: Das Array identifiziert vier separate Intervalle als innerhalb des Polygons.

Einige dieser Intervalle, z. B. das Intervall von x = 60 bis x = 63,3, schneiden keine Scheitelpunkte: Durch eine schnelle Überprüfung der x-Koordinaten aller Scheitelpunkte mit y = 10 werden solche Intervalle eliminiert.

Während des Scans können wir die Länge dieser Intervalle überwachen und dabei die Daten zu den bisher gefundenen Intervallen maximaler Länge beibehalten.

Beachten Sie einige der Auswirkungen dieses Ansatzes. Ein "v" -förmiger Scheitelpunkt ist, wenn er angetroffen wird, der Ursprung von zwei Kanten. Daher treten beim Überqueren zwei Schalter auf. Diese Schalter heben sich auf. Ein auf dem Kopf stehendes "v" wird nicht einmal verarbeitet, da beide Kanten beseitigt werden, bevor der Scan von links nach rechts gestartet wird. In beiden Fällen blockiert ein solcher Scheitelpunkt ein horizontales Segment nicht.

Mehr als zwei Kanten können einen Scheitelpunkt teilen: Dies ist bei (10,0), (60,5), (25, 20) und - obwohl es schwer zu sagen ist - bei (20,10) und (40) dargestellt 10). (Das liegt daran, dass das Baumeln geht (20,10) -> (40,10) -> (40,0) -> (40, -50) -> (40, 10) -> (20, 10) Beachten Sie, dass sich der Scheitelpunkt bei (40,0) auch im Inneren einer anderen Kante befindet ... das ist böse.) Dieser Algorithmus behandelt diese Situationen in Ordnung.

Ganz unten ist eine knifflige Situation dargestellt: Die x-Koordinaten der dortigen nichthorizontalen Segmente

30, 50

Dies bewirkt, dass alles links von x = 30 als außen betrachtet wird, alles zwischen 30 und 50 als innen und alles nach 50 wieder als außen. Der Scheitelpunkt bei x = 40 wird in diesem Algorithmus nicht einmal berücksichtigt.

So sieht das Polygon am Ende des Scans aus. Ich zeige alle vertexhaltigen Innenintervalle in Dunkelgrau, alle Intervalle mit maximaler Länge in Rot und färbe die Vertices entsprechend ihrer y-Koordinaten. Das maximale Intervall beträgt 64 Einheiten.

Die einzigen geometrischen Berechnungen bestehen darin, zu berechnen, wo Kanten horizontale Linien schneiden: das ist eine einfache lineare Interpolation. Berechnungen sind auch erforderlich, um zu bestimmen, welche inneren Segmente Scheitelpunkte enthalten: Dies sind Zwischengleichheitsbestimmungen , die leicht mit ein paar Ungleichungen berechnet werden können. Diese Einfachheit macht den Algorithmus robust und sowohl für Ganzzahl- als auch für Gleitkommakoordinatendarstellungen geeignet.

Wenn die Koordinaten geografisch sind , liegen die horizontalen Linien tatsächlich auf den Breitengraden. Ihre Längen sind nicht schwer zu berechnen: Multiplizieren Sie einfach ihre euklidischen Längen mit dem Kosinus ihres Breitengrades (in einem sphärischen Modell). Daher passt sich dieser Algorithmus gut an geografische Koordinaten an. (Um das Umwickeln der + -180-Meridianbohrung zu handhaben, muss möglicherweise zuerst eine Kurve vom Südpol zum Nordpol gefunden werden, die nicht durch das Polygon verläuft. Nachdem alle x-Koordinaten als horizontale Verschiebungen relativ dazu ausgedrückt wurden Kurve wird dieser Algorithmus das maximale horizontale Segment korrekt finden.)

Der folgende RCode wurde implementiert, um die Berechnungen durchzuführen und die Illustrationen zu erstellen.

#
# Plotting functions.
#
points.polygon <- function(p, ...) {
  points(p$v, ...)
}
plot.polygon <- function(p, ...) {
  apply(p$e, 1, function(e) lines(matrix(e[c("x.min", "x.max", "y.min", "y.max")], ncol=2), ...))
}
expand <- function(bb, e=1) {
  a <- matrix(c(e, 0, 0, e), ncol=2)
  origin <- apply(bb, 2, mean)
  delta <-  origin %*% a - origin
  t(apply(bb %*% a, 1, function(x) x - delta))
}
#
# Convert polygon to a better data structure.
#
# A polygon class has three attributes:
#   v is an array of vertex coordinates "x" and "y" sorted by increasing y;
#   e is an array of edges from (x.min, y.min) to (x.max, y.max) with y.max >= y.min, sorted by y.min;
#   bb is its rectangular extent (x0,y0), (x1,y1).
#
as.polygon <- function(p) {
  #
  # p is a list of linestrings, each represented as a sequence of 2-vectors 
  # with coordinates in columns "x" and "y". 
  #
  f <- function(p) {
    g <- function(i) {
      v <- p[(i-1):i, ]
      v[order(v[, "y"]), ]
    }
    sapply(2:nrow(p), g)
  }
  vertices <- do.call(rbind, p)
  edges <- t(do.call(cbind, lapply(p, f)))
  colnames(edges) <- c("x.min", "x.max", "y.min", "y.max")
  #
  # Sort by y.min.
  #
  vertices <- vertices[order(vertices[, "y"]), ]
  vertices <- vertices[!duplicated(vertices), ]
  edges <- edges[order(edges[, "y.min"]), ]

  # Maintaining an extent is useful.
  bb <- apply(vertices <- vertices[, c("x","y")], 2, function(z) c(min(z), max(z)))

  # Package the output.
  l <- list(v=vertices, e=edges, bb=bb); class(l) <- "polygon"
  l
}
#
# Compute the maximal horizontal interior segments of a polygon.
#
fetch.x <- function(p) {
  #
  # Update moves the line from the previous level to a new, higher level, changing the
  # state to represent all edges originating or strictly passing through level `y`.
  #
  update <- function(y) {
    if (y > state$level) {
      state$level <<- y
      #
      # Remove edges below the new level from state$current.
      #
      current <- state$current
      current <- current[current[, "y.max"] > y, ]
      #
      # Adjoin edges at this level.
      #
      i <- state$i
      while (i <= nrow(p$e) && p$e[i, "y.min"] <= y) {
        current <- rbind(current, p$e[i, ])
        i <- i+1
      }
      state$i <<- i
      #
      # Sort the current edges by x-coordinate.
      #
      x.coord <- function(e, y) {
        if (e["y.max"] > e["y.min"]) {
          ((y - e["y.min"]) * e["x.max"] + (e["y.max"] - y) * e["x.min"]) / (e["y.max"] - e["y.min"])
        } else {
          min(e["x.min"], e["x.max"])
        }
      }
      if (length(current) > 0) {
        x.array <- apply(current, 1, function(e) x.coord(e, y))
        i.x <- order(x.array)
        current <- current[i.x, ]
        x.array <- x.array[i.x]     
        #
        # Scan and mark each interval as interior or exterior.
        #
        status <- FALSE
        interior <- numeric(length(x.array))
        for (i in 1:length(x.array)) {
          if (current[i, "y.max"] == y) {
            interior[i] <- TRUE
          } else {
            status <- !status
            interior[i] <- status
          }
        }
        #
        # Simplify the data structure by retaining the last value of `interior`
        # within each group of common values of `x.array`.
        #
        interior <- sapply(split(interior, x.array), function(i) rev(i)[1])
        x.array <- sapply(split(x.array, x.array), function(i) i[1])

        print(y)
        print(current)
        print(rbind(x.array, interior))


        markers <- c(1, diff(interior))
        intervals <- x.array[markers != 0]
        #
        # Break into a list structure.
        #
        if (length(intervals) > 1) {
          if (length(intervals) %% 2 == 1) 
            intervals <- intervals[-length(intervals)]
          blocks <- 1:length(intervals) - 1
          blocks <- blocks - (blocks %% 2)
          intervals <- split(intervals, blocks)  
        } else {
          intervals <- list()
        }
      } else {
        intervals <- list()
      }
      #
      # Update the state.
      #
      state$current <<- current
    }
    list(y=y, x=intervals)
  } # Update()

  process <- function(intervals, x, y) {
    # intervals is a list of 2-vectors. Each represents the endpoints of
    # an interior interval of a polygon.
    # x is an array of x-coordinates of vertices.
    #
    # Retains only the intervals containing at least one vertex.
    between <- function(i) {
      1 == max(mapply(function(a,b) a && b, i[1] <= x, x <= i[2]))
    }
    is.good <- lapply(intervals$x, between)
    list(y=y, x=intervals$x[unlist(is.good)])
    #intervals
  }
  #
  # Group the vertices by common y-coordinate.
  #
  vertices.x <- split(p$v[, "x"], p$v[, "y"])
  vertices.y <- lapply(split(p$v[, "y"], p$v[, "y"]), max)
  #
  # The "state" is a collection of segments and an index into edges.
  # It will updated during the vertical line sweep.
  #
  state <- list(level=-Inf, current=c(), i=1, x=c(), interior=c())
  #
  # Sweep vertically from bottom to top, processing the intersection
  # as we go.
  #
  mapply(function(x,y) process(update(y), x, y), vertices.x, vertices.y)
}


scale <- 10
p.raw = list(scale * cbind(x=c(0:10,7,6,0), y=c(3,0,0,-1,-1,-1,0,-0.5,0.75,1,4,1.5,0.5,3)),
             scale *cbind(x=c(1,1,2.4,2,4,4,4,4,2,1), y=c(0,1,2,1,1,0,-0.5,1,1,0)),
             scale *cbind(x=c(6,7,6,6), y=c(.5,2,3,.5)))

#p.raw = list(cbind(x=c(0,2,1,1/2,0), y=c(0,0,2,1,0)))
#p.raw = list(cbind(x=c(0, 35, 100, 65, 0), y=c(0, 50, 100, 50, 0)))

p <- as.polygon(p.raw)

results <- fetch.x(p)
#
# Find the longest.
#
dx <- matrix(unlist(results["x", ]), nrow=2)
length.max <- max(dx[2,] - dx[1,])
#
# Draw pictures.
#
segment.plot <- function(s, length.max, colors,  ...) {
  lapply(s$x, function(x) {
      col <- ifelse (diff(x) >= length.max, colors[1], colors[2])
      lines(x, rep(s$y,2), col=col, ...)
    })
}
gray <- "#f0f0f0"
grayer <- "#d0d0d0"
plot(expand(p$bb, 1.1), type="n", xlab="x", ylab="y", main="After the Scan")
sapply(1:length(p.raw), function(i) polygon(p.raw[[i]], col=c(gray, "White", grayer)[i]))
apply(results, 2, function(s) segment.plot(s, length.max, colors=c("Red", "#b8b8a8"), lwd=4))
plot(p, col="Black", lty=3)
points(p, pch=19, col=round(2 + 2*p$v[, "y"]/scale, 0))
points(p, cex=1.25)

Hier ist eine rasterbasierte Lösung. Es ist schnell (ich habe die ganze Arbeit von Anfang bis Ende in 14 Minuten erledigt), erfordert kein Scripting, nimmt nur wenige Operationen und ist einigermaßen genau.

Beginnen Sie mit einer Rasterdarstellung des Polygons. Dieser verwendet ein Raster von 550 Zeilen und 1200 Spalten:

In dieser Darstellung haben die grauen (inneren) Zellen den Wert 1 und alle anderen Zellen sind NoData.

Berechnen Sie die Durchflussakkumulation in West-Ost-Richtung unter Verwendung der Einheitszellenwerte für das Gewichtsgitter (Menge des "Niederschlags"):

Niedrige Akkumulation ist dunkel und steigt zu höchsten Akkumulationen im hellen Gelb an.

Ein zonales Maximum (unter Verwendung des Polygons für das Gitter und der Flussakkumulation für die Werte) kennzeichnet die Zelle (n), in der / denen der Fluss am größten ist. Um dies zu zeigen, musste ich nach rechts unten zoomen:

Die roten Zellen markieren die Enden der höchsten Fließansammlungen: Sie sind die äußersten rechten Endpunkte der Innensegmente maximaler Länge des Polygons.

Um diese Segmente zu finden, platzieren Sie das gesamte Gewicht auf den roten Blutkörperchen und leiten Sie den Fluss rückwärts!

Der rote Streifen am unteren Rand markiert zwei Zellenreihen: Darin befindet sich das horizontale Segment maximaler Länge. Verwenden Sie diese Darstellung unverändert für die weitere Analyse oder konvertieren Sie sie in eine Polylinien- (oder Polygon-) Form.

Bei einer Raster-Darstellung ist ein Diskretisierungsfehler aufgetreten. Sie kann durch Erhöhen der Auflösung mit einem gewissen Aufwand an Rechenzeit verringert werden.

Ein wirklich schöner Aspekt dieses Ansatzes ist, dass wir in der Regel extreme Werte von Dingen als Teil eines größeren Workflows finden, in dem ein Ziel erreicht werden muss: das Aufstellen einer Pipeline oder eines Fußballfelds, das Erstellen von ökologischen Puffern und so weiter. Der Prozess beinhaltet Kompromisse. Daher ist die längste horizontale Linie möglicherweise nicht Teil einer optimalen Lösung. Wir könnten uns stattdessen darum kümmern, zu wissen, wo fast die längsten Zeilen liegen würden. Dies ist einfach: Anstatt den zonalen Maximalfluss auszuwählen, markieren Sie alle Zellen, die sich in der Nähe eines zonalen Maximums befinden. In diesem Beispiel entspricht das zonale Maximum 744 (die Anzahl der Spalten, die vom längsten inneren Segment überspannt werden). Wählen Sie stattdessen alle Zellen innerhalb von 5% dieses Maximums aus:

Läuft der Fluss von Ost nach West, entsteht diese Sammlung horizontaler Segmente:

Dies ist eine Karte von Orten, an denen die ununterbrochene Ost-West-Ausdehnung mindestens 95% der maximalen Ost-West-Ausdehnung innerhalb des Polygons beträgt.

In Ordnung. Ich habe eine andere (bessere) Idee ( Idee-Nr. 2 ). Aber ich nehme an, dass es besser ist, als Python-Skript, nicht als SQL-Querry, realisiert zu werden. Auch hier ist der häufigste Fall, nicht nur EW.

Sie benötigen einen Begrenzungsrahmen für das Polygon und einen Azimut (A) als Messrichtung. Angenommen, die Länge der BBox-Kanten ist LA und LB. Der maximal mögliche Abstand (MD) innerhalb eines Polygons beträgt:MB = (LA^2 * LB^2)^(1/2) , so Wert suchen (V) ist nicht größer als MB: V <= MB.

1. Erstellen Sie ausgehend von einem beliebigen Scheitelpunkt der BBox eine Linie (LL) mit der Länge MB und dem Azimut A.
2. Schneide die Linie LL mit dem Polygon, um die Schnittlinie (IL) zu erhalten.
3. Überprüfen Sie die Geometrie von IL - wenn es nur zwei Punkte in der IL-Linie gibt, berechnen Sie die Länge. Wenn 4 oder mehr - berechnen Sie die Segmente und ermitteln Sie die Länge des längsten. Null (überhaupt keine Kreuzung) - ignorieren.
4. Erstellen Sie weitere LL-Linien, die sich vom Startpunkt aus gegen den Uhrzeigersinn in Richtung der Kanten der BBox bewegen, bis Sie nicht mehr am Startpunkt enden (Sie machen die gesamte Schleife über die BBox).
5. Nehmen Sie den größten IL-Längenwert auf (tatsächlich müssen Sie nicht alle Längen speichern, Sie können beim Schleifen nur den "bisherigen" Maximalwert beibehalten) - es wird das sein, wonach Sie suchen.

Ich bin mir nicht sicher, ob Fetzers Antwort das ist, was Sie tun möchten, aber es ist so, dass die st_box2d den Job erledigen kann.

Die Idee von SS_Rebelious Nr. 1 wird in vielen Fällen funktionieren, jedoch nicht für einige konkave Polygone.

Ich denke, dass Sie künstliche LW-Linien erstellen müssen, deren Punkte Kanten folgen, wenn Scheitelpunktlinien die Grenzen des Polygons überschreiten, wenn eine Ost-West-Linienmöglichkeit besteht.

Dazu können Sie versuchen, ein Polygon mit 4 Knoten zu erstellen, bei dem die Linienlänge hoch ist, und das Polygon P * erstellen, bei dem es sich um das vorhergehende Polygon handelt, und prüfen, ob die Min (y1) und Max (y2) eine X-Linie hinterlassen Möglichkeit. (wobei y1 die Menge der Punkte zwischen der linken oberen Ecke und der rechten oberen Ecke ist und y2 die Menge von y zwischen der linken unteren Ecke und der rechten unteren Ecke Ihres 4-Knoten-Polygons). Das ist nicht so einfach. Ich hoffe, Sie finden psql-Tools, die Ihnen helfen!

Ich habe eine Idee-№1 ( Bearbeiten: für den allgemeinen Fall, nicht nur EW Richtung, und mit einigen Einschränkungen, die in den Kommentaren beschrieben sind). Ich werde den Code nicht zur Verfügung stellen, nur ein Konzept. Die "x-Richtung" ist tatsächlich ein Azimut, der durch ST_Azimut berechnet wird. Die vorgeschlagenen Schritte sind:

1. Extrahieren Sie alle Eckpunkte aus dem Polygon als Punkte.
2. Erstellen Sie Linien zwischen jedem Punktepaar.
3. Wählen Sie Linien aus (nennen wir sie LW-Linien), die sich innerhalb des ursprünglichen Polygons befinden (wir benötigen keine Linien, die die Grenzen des Polygons überschreiten).
4. Finde Entfernungen und Azimute für jede LW-Linie.
5. Wählen Sie den größten Abstand von LW-Linien, bei denen Azimut gleich gesuchtem Azimut ist oder in einem bestimmten Intervall liegt (es kann sein, dass kein Azimut genau gleich gesuchtem Azimut ist).

Schauen Sie sich meine Frage und die Antwort an von Evil Genius an.

Hoffentlich hat Ihr See-Polygon mehrere Punkte. Sie können auf diesen Punkten Linien mit einem Azimut (Aspekt, Georichtung) erstellen. Wählen Sie die Länge der ausreichend großen Linien (ST_MakePoint-Teil), damit Sie die kürzeste Linie zwischen den beiden am weitesten entfernten Linien berechnen können.

Hier ist ein Beispiel:

Das Beispiel zeigt die maximale Breite des Polygons. Für diesen Ansatz wähle ich ST_ShortestLine (rote Linie). ST_MakeLine würde den Wert erhöhen (blaue Linie) und der Endpunkt der Linie (unten links) würde die blaue Linie des Polygons treffen. Sie müssen den Abstand mit den Schwerpunkten der erstellten (Hilfe-) Linien berechnen.

Eine Idee für unregelmäßige oder konkave Polygone für diesen Ansatz. Möglicherweise müssen Sie das Polygon mit einem Raster schneiden.