PostGIS-Polygonkantenanalyse (Ausrichtung, Kantenlänge)

9

Ich bin ziemlich neu in der Welt von GIS und insbesondere von PostGIS. Entschuldigen Sie mich bitte, wenn die Antwort offensichtlich erscheint ...

Ich möchte einige Gebäude analysieren. Eine Sache, die mich interessiert, sind ihre Fassadenoberflächen zusammen mit der jeweiligen Ausrichtung. Wie im Bild unten dargestellt, möchte ich die Länge und (normale) Ausrichtung aller Kanten in einer Reihe von Polygonen haben. Im Beispiel habe ich nur eine Oberfläche hervorgehoben.

Geben Sie hier die Bildbeschreibung ein

Eine Ergebnistabelle könnte folgendermaßen aussehen:

building_id | edge_id | orientation | edge_length
-------------------------------------------------
      1     |    1    |     315     |    10.0
      1     |    2    |      45     |     7.0
      1     |   ...   |     ...     |     ...

Ich bin mir jedoch nicht sicher, ob es eine clevere Möglichkeit ist, das Ergebnis für die weitere Verarbeitung zu speichern (z. B. Abstand von der Kante zum nächsten Gebäude berechnen usw.). Meine Frage ist also zweifach:

  1. Gibt es eine effiziente PostGIS-Funktion, mit der die Kanten eines Polygons analysiert werden können? Falls es keine native PostGIS-Funktion gibt, wäre ich alternativ an einem Python-basierten Ansatz interessiert.
  2. Was wäre eine clevere Möglichkeit, das Ergebnis in einer PostGIS-Tabelle zu speichern, da die Polygone möglicherweise eine unterschiedliche Anzahl von Kanten haben?
n1000
quelle
2
Erstellen Sie zuerst die Segmente des Polygons: stackoverflow.com/questions/7595635/… Dann sollten der Startpunkt und die Endpunktkoordinaten zu Spalten wie x1, y1 und x2, y2 und dann zu ST_Azimuth (ST_Startpoint (Geometrien), ST_Endpoint (Geometrien)) gehen. . ( postgis.org/docs/ST_Azimuth.html )
Tamas Kosa
1
@TamasKosa: Du hast die Essenz einer guten Antwort. Warum nicht zu einem erweitern? Für Normalen benötigen Azimute +/- pi / 2.
Martin F
1
@TamasKosa Dies ist ein Ansatz, über den ich auch nachgedacht habe. Verwenden Sie ST_ExteriorRing und erhalten Sie dann die Azimute, wie Sie sagen. Wie würde ich die Ergebnisse ideal speichern, da Gebäude unterschiedlich viele Kanten haben können? In einer Tabelle wie oben beschrieben? Ich stimme MartinF zu, das ist fast eine Antwort;)
n1000
Nur neugierig, warum willst du Normalen ... Sonneneinstrahlung?
Martin F
1
Der erste Teil Ihrer Frage wurde beantwortet - Sie können ST_Dumppoints und ST_Azimuth verwenden. Für den zweiten Teil, da die Ausgabe keine räumlichen Elemente enthält, würde ich einen Link zu polygonID und edge_id denken, wie Sie ihn finden würden.
John Powell

Antworten:

10

Gestern hatte ich keine Zeit, es im Detail zu erstellen ... Sehen Sie meine Lösung in 4 Schritten:

CREATE OR REPLACE VIEW bd_segment AS
SELECT
      ST_PointN(geom, generate_series(1, ST_NPoints(geom)-1)) AS sp,
      ST_PointN(geom, generate_series(2, ST_NPoints(geom)  )) AS ep
    FROM
       -- extract the individual linestrings
       (SELECT (ST_Dump(ST_Boundary(the_geom))).geom
       FROM bd) AS linestrings;

CREATE OR REPLACE VIEW bd_segment_geom AS
SELECT sp, ep, st_makeline(sp, ep) 
FROM bd_segment;

CREATE OR REPLACE view bd_segment_id AS 
SELECT bd.gid, row_number() 
    OVER (order by bd.gid), degrees(st_azimuth(ff.sp, ff.ep)-1.57079633) AS az_deg,
    ST_LENGTH(ff.st_makeline) , ff.st_makeline FROM bd_segment_geom ff
JOIN bd ON st_touches(ff.st_makeline, bd.the_geom)
GROUP BY bd.gid, ff.sp, ff.ep, ff.st_makeline;

UPDATE bd_segment_id
SET az_deg = az_deg + 360
WHERE az_deg < 0;

Bei der letzten Abfrage erhalten Sie die Gebäude-IDs mit einer räumlichen Verknüpfung mithilfe von st_touches. Ich hoffe es hilft. Update - In qgis sieht die Lösung folgendermaßen aus: Geben Sie hier die Bildbeschreibung ein

Tamas Kosa
quelle
1
Beeindruckend! Ich habe es zum Laufen gebracht. Ich danke dir sehr! Bei einer großen Anzahl von Gebäuden wird es etwas langsam. Die Azimute sind momentan keine Normalen. Hätten Sie auch eine Idee, wie Sie das angehen können? Ich bin nicht sicher, wie ich die "äußere" Seite des Polygons finden soll.
n1000
1
Addiere 90 Grad zum Azimut im Bogenmaß wie folgt: Grad (st_azimuth (ff.sp, ff.ep) +1.57079633). Manchmal werden Werte generiert, die größer als 360 sind. Mit einer Aktualisierungsabfrage können Sie diese jedoch ersetzen. Wenn Sie es als statische Ansicht verwenden möchten, erstellen Sie "MATERIALISIERTE ANSICHT ERSTELLEN" und es wird nicht gleich beim ersten Mal langsam.
Tamas Kosa
2
Nicht ganz. Angenommen, Nord ist 0 °, ergibt dies die Normalen zur Innenseite des Polygons / Gebäudes (wie auch in Ihrem Screenshot zu sehen). Aber Sie haben Recht - ein einfacher UPDATEsollte den Trick tun. Nochmals vielen Dank für diese großartige Lösung. Ich werde noch einige Tage warten, wenn andere Antworten erscheinen, bevor ich akzeptiere.
n1000
1
Was ist mit ST_ForceRHR? Diese Antwort scheint tatsächlich in Ordnung zu sein.
Jakub Kania
@JakubKania Ich habe versucht, eine ST_ForceRHRLösung zu finden, war aber nicht erfolgreich. Wäre dankbar für Hinweise ... Ich habe versuchtST_Dump(ST_Boundary(ST_ForceRHR(the_geom)))
n1000