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

13

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.

Kardiopteris
quelle
3
Ich bin mit der Postgis-Funktionalität nicht vertraut. Möglicherweise ist jedoch ein Begrenzungsrahmen-Werkzeug vorhanden. Die Breite des Begrenzungsrahmens ist der maximale Abstand in EW-Richtung.
Fezter
4
@Fetzter das stimmt nicht: ein Gegenbeispiel, auch für ein einfaches komplexes Polygon, ist eine dünne Raute, die sich von SW nach NE erstreckt. Seine maximale Ost-West-Breite kann ein beliebig kleiner Bruchteil der Breite seines Begrenzungsrahmens sein.
Whuber
1
Ich habe für diese Aufgabe ein Hilfsprogramm erstellt, das auf diesen und diesen Vorschlägen basiert . Es kann die maximale oder minimale Breite des Polygons berechnet werden. Momentan funktioniert es mit shp-Dateien, aber Sie können es umschreiben, um mit einem PostGIS zu arbeiten, oder einfach eine Weile warten, bis es sich zu einem QGIS-Plugin entwickelt, das auch mit PostGIS funktioniert. Detaillierte Beschreibung und Download-Link hier .
SS_Rebelious

Antworten:

16

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

Ein Polygon

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.

Nach dem Scan

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)
whuber
quelle
Gibt es einen Satz, der beweist, dass die Linie maximaler Länge innerhalb eines nicht konvexen Polygons in einer bestimmten Richtung mindestens einen Scheitelpunkt dieses Polygons schneidet?
SS_Rebelious
@SS Ja, das gibt es. Hier ist eine Beweisskizze: Wenn es keine Schnittmenge gibt, liegen die Endpunkte des Segments beide auf den Innenkanten, und das Segment kann zumindest ein wenig nach oben und unten bewegt werden. Seine Länge ist eine lineare Funktion des Verschiebungsbetrags. Daher kann es nur dann eine maximale Länge haben, wenn sich die Länge beim Verschieben nicht ändert. Dies impliziert sowohl (a) als auch, dass es Teil eines Parallelogramms ist, das aus Segmenten maximaler Länge besteht, und (b) dass sowohl die Ober- als auch die Unterkante dieses Parallelogramms einen Scheitelpunkt QED erfüllen müssen.
Whuber
Und wie heißt dieser Satz? Ich kämpfe darum, es zu finden. Übrigens, was ist mit gekrümmten Kanten ohne Scheitelpunkt (ich meine einen theoretischen Ansatz)? Eine Skizze eines Beispiels der Figur, die ich meine (ein hantelförmiges Polygon): "C = D".
SS_Rebelious
@SS Wenn die Kanten gekrümmt sind, gilt der Satz nicht mehr. Techniken der Differentialgeometrie können angewendet werden, um nützliche Ergebnisse zu erhalten. Ich habe diese Methoden aus dem Buch von Cheeger & Ebin, Vergleichssätze in der Riemannschen Geometrie, gelernt . Da die meisten GISs Kurven jedoch ohnehin durch detaillierte Polylinien approximieren, ist die Frage (aus praktischen Gründen) umstritten.
whuber
Könnten Sie den Namen des Theorems (und die Seite, wenn möglich) angeben? Ich habe das Buch und konnte den benötigten Satz nicht finden.
SS_Rebelious
9

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:

Polygon

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

Fließansammlung

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:

Maximal

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!

Ergebnis

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:

Ausgewählte nahezu optimale Zellen

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

Nahezu optimale Lösungen

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.

whuber
quelle
3

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.
SS_Rebelious
quelle
Dies klingt wie eine Doppelschleife über den Eckpunkten: Dies ist ineffizient genug, um vermieden zu werden (mit Ausnahme stark vereinfachter Polygone).
Whuber
@whuber, ich sehe hier keine zusätzlichen Schleifen. Es gibt nur einige bedeutungslose Verarbeitungen von 2 Seiten von BB, die nur Nullen ergeben. Diese Verarbeitung wurde jedoch in dem Skript, das ich in der hier gelöschten Antwort angegeben habe, ausgeschlossen (es scheint, dass es sich derzeit um einen Kommentar handelt, ich sehe ihn jedoch nicht als Kommentar - nur als gelöschte Antwort).
SS_Rebelious
(1) Dies ist der dritte Kommentar zur Frage. (2) Sie haben Recht: Wenn Sie Ihre Beschreibung sehr sorgfältig lesen, erscheint es mir jetzt, als ob Sie das längste Segment zwischen (den vier) Eckpunkten des Begrenzungsrahmens und den Eckpunkten des Polygons finden. Ich verstehe jedoch nicht, wie dies die Frage beantwortet: Das Ergebnis ist definitiv nicht das, wonach das OP gesucht hat.
Whuber
@whuber, vorgeschlagener Algorithmus findet den längsten Schnittpunkt des Polygons mit der Linie, die die angegebene Richtung darstellt. Anscheinend IST das Ergebnis, was gefragt wurde, ob der Abstand zwischen Schnittlinien -> 0 ist oder alle Eckpunkte (für nicht gekrümmte Figuren) passiert.
SS_Rebelious
3

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. Ein Beispiel, wo es nicht funktioniert

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!

Ein Name
quelle
Das ist auf dem richtigen Weg. Das längste EW-Segment befindet sich unter den Schnittpunkten mit dem Inneren des Polygons, wobei horizontale Linien durch die Eckpunkte des Polygons verlaufen. Dies erfordert Code, um die Eckpunkte zu durchlaufen. Es gibt eine alternative (aber äquivalente) Methode, bei der eine künstliche Ost-West-Strömung über eine Rasterdarstellung des Polygons verfolgt wird: Die maximale Strömungslänge im Polygon (eine seiner "Zonenstatistiken") entspricht der gewünschten Breite. Die Rasterlösung wird in nur 3 oder 4 Schritten erstellt und erfordert keine Schleifen oder Skripte.
Whuber
@Aname, bitte fügen Sie "№1" zur "Idee von SS_Rebelious" hinzu, um Missverständnisse zu vermeiden: Ich habe einen weiteren Vorschlag hinzugefügt. Ich kann Ihre Antwort nicht selbst bearbeiten, da diese Bearbeitung weniger als 6 Zeichen umfasst.
SS_Rebelious
1

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).
SS_Rebelious
quelle
Dies funktioniert auch bei einigen Dreiecken nicht , z. B. bei den Eckpunkten (0,0), (1000, 1000) und (501, 499). Seine maximale Ost-West-Breite beträgt ungefähr 2; Die Azimute sind alle um 45 Grad; und unabhängig davon ist das kürzeste Liniensegment zwischen Eckpunkten mehr als 350-mal länger als die Ost-West-Breite.
Whuber
@whuber, Sie haben Recht, es wird für Dreiecke scheitern, aber für Polygone, die einige Naturmerkmale darstellen, kann es nützlich sein.
SS_Rebelious
1
Es ist schwer, ein Verfahren zu empfehlen, das selbst in einfachen Fällen dramatisch scheitert, in der Hoffnung, dass es manchmal eine richtige Antwort erhält!
whuber
@whuber, also empfehle es nicht! ;-) Ich habe diese Problemumgehung vorgeschlagen, da es keine Antworten auf diese Frage gab. Beachten Sie, dass Sie möglicherweise Ihre eigene bessere Antwort veröffentlichen. Übrigens, wenn Sie einige Punkte auf die Dreiecksränder setzen, wird mein Vorschlag funktionieren ;-)
SS_Rebelious
Ich habe mehrere Ansätze vorgeschlagen. Ein Raster ist bei gis.stackexchange.com/questions/32552/... und ausgearbeitet an forums.esri.com/Thread.asp?c=93&f=982&t=107703&mc=3 . Eine andere - nicht ganz zutreffend, aber mit interessanten Verwendungsmöglichkeiten - finden Sie unter gis.stackexchange.com/questions/23664/… (die Radon-Transformation). Dies ist unter stats.stackexchange.com/a/33102 dargestellt .
Whuber
1

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:

Bildbeschreibung hier eingeben

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.

Stefan
quelle