Análisis del borde del polígono PostGIS (orientación, longitud del borde)

9

Soy bastante nuevo en el mundo de los SIG y especialmente en PostGIS, así que discúlpeme si la respuesta parece evidente ...

Me gustaría hacer un análisis de varios edificios. Una cosa que me interesa es su superficie de fachada junto con la orientación respectiva. Como se ilustra en la imagen a continuación, me gustaría tener la longitud y la orientación (normal) de todos los bordes en una serie de polígonos. En el ejemplo, destaqué solo una superficie.

ingrese la descripción de la imagen aquí

Una tabla de resultados podría verse así:

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

Sin embargo, no estoy seguro de si es una forma inteligente de almacenar el resultado para su posterior procesamiento (por ejemplo, calcular la distancia desde el borde hasta el próximo edificio, etc.). Entonces mi pregunta es doble:

  1. ¿Existe una función eficiente de PostGIS que pueda analizar los bordes de un polígono? En caso de que no haya una función nativa de PostGIS, alternativamente me interesaría un enfoque basado en Python.
  2. ¿Cuál sería una forma inteligente de almacenar el resultado en una tabla PostGIS, ya que los polígonos pueden tener diferentes números de bordes?
n1000
fuente
2
Primero cree los segmentos del polígono: stackoverflow.com/questions/7595635/… Luego, el punto inicial y las coordenadas del punto final deben ir a columnas como x1, y1 y x2, y2 y que ST_Azimuth (ST_Startpoint (geometrías), ST_Endpoint (geometrías)) . ( postgis.org/docs/ST_Azimuth.html )
Tamas Kosa
1
@TamasKosa: Tienes la esencia de una buena respuesta. ¿Por qué no expandirlo en uno? Además, para las normales, los acimutes necesitan +/- pi / 2.
Martin F
1
@TamasKosa Este es un enfoque en el que también estaba pensando. Use ST_ExteriorRing y luego obtenga los acimutes como usted dice. ¿Cómo almacenaría idealmente los resultados, ya que los edificios pueden tener un número diferente de bordes? En una tabla como la que describí anteriormente? Estoy de acuerdo con MartinF, esta es casi una respuesta;)
n1000
Solo por curiosidad, ¿por qué quieres normales ... exposición al sol?
Martin F
1
La primera parte de su pregunta ha sido respondida: puede usar ST_Dumppoints y ST_Azimuth. Para la segunda parte, como no hay elementos espaciales en la salida, creo que sería encontrar un enlace a polygonID y edge_id como usted tiene.
John Powell

Respuestas:

10

Ayer no tuve tiempo para crearlo en detalles ... Vea mi solución en 4 pasos:

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 última consulta le proporciona los identificadores de construcción con una unión espacial usando st_touches. Espero eso ayude. Actualización: en qgis la solución se ve así: ingrese la descripción de la imagen aquí

Tamas Kosa
fuente
1
¡Impresionante! Yo tengo que trabajar. Muchas gracias! Se vuelve un poco lento con una gran cantidad de edificios. Los acimutes no son normales en este momento. ¿Tendrías también una idea de cómo abordar eso? No estoy seguro de cómo encontrar el lado "exterior" del polígono.
n1000
1
agregue 90 grados al azimut en radianes como este: grados (st_azimuth (ff.sp, ff.ep) +1.57079633). A veces generará valores superiores a 360. pero con una consulta de actualización puede reemplazarlos. Si desea utilizarlo como una vista estática, cree "CREAR VISTA MATERIALIZADA" y no será lento solo la primera vez.
Tamas Kosa
2
No exactamente. Asumiendo que el Norte es 0 °, esto dará la normalidad hacia el interior del polígono / edificio (como también se ve en la captura de pantalla). Pero tienes razón: un simple UPDATEdebería hacer el truco. Gracias de nuevo por esta gran solución. Esperaré algunos días más si aparecen otras respuestas, antes de aceptar.
n1000
1
¿Qué hay de ST_ForceRHR? Esta respuesta en realidad parece estar bien.
Jakub Kania
@JakubKania Intenté encontrar una ST_ForceRHRsolución, pero no tuve éxito. Estaría agradecido por las sugerencias ... Lo intentéST_Dump(ST_Boundary(ST_ForceRHR(the_geom)))
n1000