Analiza krawędzi wielokąta PostGIS (orientacja, długość krawędzi)

Jestem raczej nowy w świecie GIS, a zwłaszcza PostGIS, więc przepraszam, jeśli odpowiedź wydaje się oczywista ...

Chciałbym przeprowadzić analizę szeregu budynków. Interesuje mnie ich powierzchnia elewacji wraz z odpowiednią orientacją. Jak pokazano na poniższym obrazku, chciałbym mieć długość i (normalną) orientację wszystkich krawędzi w szeregu wielokątów. W przykładzie wyróżniłem tylko jedną powierzchnię.

wprowadź opis zdjęcia tutaj

Tabela wyników może wyglądać następująco:

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

Nie jestem jednak pewien, czy jest to inteligentny sposób na przechowywanie wyników do dalszego przetwarzania (np. Oblicz odległość od krawędzi do następnego budynku itp.). Więc moje pytanie jest dwojakie:

Czy istnieje skuteczna funkcja PostGIS, która może analizować krawędzie wielokąta? W przypadku, gdy nie ma natywnej funkcji PostGIS, alternatywnie byłbym zainteresowany podejściem opartym na języku Python.
Jaki byłby inteligentny sposób przechowywania wyniku w tabeli PostGIS, ponieważ wielokąty mogą mieć różną liczbę krawędzi?

— n1000
źródło

Najpierw utwórz segmenty wielokąta: stackoverflow.com/questions/7595635/… Następnie współrzędne punktu początkowego i końcowego powinny przejść do kolumn takich jak x1, y1 i x2, y2 i niż ST_Azimuth (ST_Startpoint (geometria), ST_Endpoint (geometria)) . ( postgis.org/docs/ST_Azimuth.html )

— Tamas Kosa

@TamasKosa: Masz istotę dobrej odpowiedzi. Dlaczego nie rozwinąć go w jedno? Ponadto, dla normalnych, azymuty potrzebują +/- pi / 2.

— Martin F

@TamasKosa O tym podejściu również myślałem. Użyj ST_ExteriorRing, a następnie uzyskaj azymuty, jak mówisz. Jak najlepiej przechowywać wyniki, ponieważ budynki mogą mieć różną liczbę krawędzi? W tabeli jak opisałem powyżej? Zgadzam się z MartinFem, to prawie odpowiedź;)

— n1000

Ciekawe, dlaczego chcesz normalnych ... ekspozycji na słońce?

— Martin F

Odpowiedź na pierwszą część pytania - możesz użyć ST_Dumppoints i ST_Azimuth. W drugiej części, ponieważ nie ma żadnych elementów przestrzennych na wyjściu, pomyślałbym, że można by znaleźć link do polygonID i edge_id.

— John Powell,

Wczoraj nie miałem czasu na szczegółowe tworzenie ... Zobacz moje rozwiązanie w 4 krokach:

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;

Ostatnie zapytanie daje identyfikatory budynku z łączeniem przestrzennym za pomocą st_touches. Mam nadzieję, że to pomoże. Aktualizacja - w qgis rozwiązanie wygląda następująco: wprowadź opis zdjęcia tutaj

— Tamas Kosa
źródło

Imponujący! Mam to do pracy. Dziękuję bardzo! Staje się trochę powolny przy dużej liczbie budynków. Azymuty nie są teraz normalne. Czy masz również pomysł, jak to rozwiązać? Nie jestem pewien, jak znaleźć „zewnętrzną” stronę wielokąta.

— n1000

dodaj 90 stopni do azymutu w radianach w następujący sposób: stopnie (st_azimuth (ff.sp, ff.ep) +1.57079633). Będzie generować czasami wartości większe niż 360. ale za pomocą zapytania o aktualizację możesz je zastąpić. Jeśli chcesz użyć jako widoku statycznego, utwórz „UTWÓRZ MATERIALIZOWANY WIDOK” i nie będzie on wolny za pierwszym razem.

— Tamas Kosa

Nie do końca. Zakładając, że północ wynosi 0 °, da to normalność w kierunku wnętrza wielokąta / budynku (jak pokazano również na zrzucie ekranu). Ale masz rację - proste UPDATEpowinno wystarczyć . Jeszcze raz dziękuję za to świetne rozwiązanie. Poczekam jeszcze kilka dni, jeśli pojawią się inne odpowiedzi, zanim zaakceptuję.

— n1000

Co ST_ForceRHR? Ta odpowiedź wydaje się w porządku.

— Jakub Kania

@JakubKania Próbowałem znaleźć ST_ForceRHRrozwiązanie, ale nie udało się. Byłbym wdzięczny za podpowiedzi ... PróbowałemST_Dump(ST_Boundary(ST_ForceRHR(the_geom)))

— n1000