Analyse de bord de polygone PostGIS (orientation, longueur de bord)

9

Je suis plutôt nouveau dans le monde des SIG et surtout PostGIS, alors veuillez m'excuser si la réponse semble évidente ...

Je voudrais faire une analyse sur un certain nombre de bâtiments. Une chose qui m'intéresse est leurs surfaces de façade ainsi que l'orientation respective. Comme illustré dans l'image ci-dessous, je voudrais avoir la longueur et l'orientation (normale) de toutes les arêtes dans une série de polygones. Dans l'exemple, j'ai mis en évidence une seule surface.

entrez la description de l'image ici

Un tableau de résultats pourrait ressembler à ceci:

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

Cependant, je ne sais pas s'il s'agit d'un moyen intelligent de stocker le résultat pour un traitement ultérieur (par exemple, calculer la distance entre le bord et le bâtiment suivant, etc.). Ma question est donc double:

  1. Existe-t-il une fonction PostGIS efficace qui peut analyser les bords d'un polygone? Dans le cas où il n'y a pas de fonction PostGIS native, je serais également intéressé par une approche basée sur Python.
  2. Quelle serait une façon intelligente de stocker le résultat dans une table PostGIS, car les polygones peuvent avoir différents nombres d'arêtes?
n1000
la source
2
Créez d'abord les segments du polygone: stackoverflow.com/questions/7595635/… Ensuite, le point de départ et les coordonnées du point final doivent aller dans des colonnes telles que x1, y1 et x2, y2 et que ST_Azimuth (ST_Startpoint (géométries), ST_Endpoint (géométries)) . ( postgis.org/docs/ST_Azimuth.html )
Tamas Kosa
1
@TamasKosa: Vous avez l'essence d'une bonne réponse. Pourquoi ne pas l'étendre en un seul? De plus, pour les normales, les azimuts ont besoin de +/- pi / 2.
Martin F
1
@TamasKosa C'est une approche à laquelle je pensais également. Utilisez ST_ExteriorRing , puis obtenez les azimuts comme vous le dites. Comment enregistrer idéalement les résultats, car les bâtiments peuvent avoir un nombre différent d'arêtes? Dans un tableau comme je l'ai décrit ci-dessus? Je suis d'accord avec MartinF, c'est presque une réponse;)
n1000
Juste curieux, pourquoi voulez-vous des expositions normales ... au soleil?
Martin F
1
La première partie de votre question a été répondue - vous pouvez utiliser ST_Dumppoints et ST_Azimuth. Pour la deuxième partie, comme il n'y a pas d'éléments spatiaux dans la sortie, je pense qu'un lien vers polygonID et edge_id comme vous le feriez serait trouvé.
John Powell

Réponses:

10

Hier je n'ai pas eu le temps de la créer en détails ... Voir ma solution en 4 étapes:

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;

La dernière requête vous donne les identifiants de construction avec une jointure spatiale à l'aide de st_touches. J'espère que cela aide. Mise à jour - Dans qgis, la solution ressemble à ceci: entrez la description de l'image ici

Tamas Kosa
la source
1
Impressionnant! Je l'ai fait fonctionner. Merci beaucoup! Cela devient un peu lent avec un grand nombre de bâtiments. Les azimuts ne sont pas normaux en ce moment. Auriez-vous également une idée de comment y remédier? Je ne sais pas comment trouver le côté "extérieur" du polygone.
n1000
1
ajoutez 90 degrés à l'azimut en radian comme ceci: degrés (st_azimuth (ff.sp, ff.ep) +1,57079633). Il générera parfois des valeurs supérieures à 360. mais avec une requête de mise à jour, vous pouvez les remplacer. Si vous souhaitez utiliser une vue statique, créez "CRÉER UNE VUE MATÉRIALISÉE" et cela ne sera pas lent juste à la première fois.
Tamas Kosa
2
Pas assez. En supposant que le nord est de 0 °, cela donnera la normale vers l' intérieur du polygone / bâtiment (comme on le voit également dans votre capture d'écran). Mais vous avez raison - un simple UPDATEdevrait faire l'affaire. Merci encore pour cette excellente solution. J'attendrai encore quelques jours si d'autres réponses apparaissent, avant d'accepter.
n1000
1
Et alors ST_ForceRHR? Cette réponse semble en fait correcte.
Jakub Kania
@JakubKania J'ai essayé de trouver une ST_ForceRHRsolution, mais sans succès. Serait reconnaissant pour des indices ... J'ai essayéST_Dump(ST_Boundary(ST_ForceRHR(the_geom)))
n1000