# PostGIS polygon edge analysis (orientation, edge length)

I am rather new to the world of GIS and especially PostGIS, so please excuse me if the answer seems evident...

I would like to do analysis on a number of buildings. One thing I am interested in is their facade surfaces along with the respective orientation. As illustrated in the picture below, I would like to have the length and (normal) orientation of all edges in a series of polygons. In the example I highlighted only one surface.

A result table could look like this:

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

However, I am not sure if it is a smart way to store the result for further processing (e.g. calculate distance from edge to next building, etc.). So my question is twofold:

1. Is there an efficient PostGIS function that can analyze a polygon's edges? In case there is no native PostGIS function I would alternatively be interested in an Python-based approach.
2. What would be a smart way to store the result in a PostGIS table, since the polygons may have different numbers of edges?
• First create the segments of the polygon: stackoverflow.com/questions/7595635/… Then the startpoint and the endpoint coordinates should go to columns like x1,y1 and x2,y2 and than ST_Azimuth(ST_Startpoint(geometries), ST_Endpoint(geometries)). (postgis.org/docs/ST_Azimuth.html) May 18, 2015 at 15:34
• @TamasKosa: You have the essence of a good answer. Why not expand it into one? Also, for normals, azimuths need +/- pi/2. May 18, 2015 at 17:33
• @TamasKosa This is an approach I was also thinking about. Use ST_ExteriorRing and then get the azimuths as you say. How would I ideally store the results, since buildings can have different number of edges? In a table like I described above? I agree with MartinF, this is almost an answer ;) May 18, 2015 at 19:03
• Just curious, why do you want normals ... sun exposure? May 19, 2015 at 0:46
• The first part of your question has been answered -- you can use ST_Dumppoints and ST_Azimuth. For the second part, as there are no spatial elements to the output, I would think a link to polygonID and edge_id as you have would be find. May 19, 2015 at 6:32

Yesterday I had no time to create it in details... See my solution in 4 steps:

``````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;
``````

The last query give you the building ids with a spatial join using st_touches. Hope it helps. Update - In qgis the solution looks like this:

• Impressive! I got it to work. Thank you so much! It becomes a little slow with a large number of buildings. The azimuths are not normals right now. Would you also have an idea how to address that? I am not sure how to find the "exterior" side of the polygon. May 19, 2015 at 10:19
• add 90 deggrees to the azimuth in radian like this: degrees(st_azimuth(ff.sp, ff.ep)+1.57079633). It will generate sometimes values that bigger than 360. but with an update query you can replace those. If you would like to use as a static view create "CREATE MATERIALIZED VIEW" and it won't be slow just at the first time. May 19, 2015 at 11:01
• Not quite. Assuming North is 0°, this will give the normal towards the inside of the polygon / building (as also seen in your screenshot). But you are right - a simple `UPDATE` should do the trick. Thanks again for this great solution. I will wait some more days if other answers appear, before accepting. May 19, 2015 at 11:46
• What about `ST_ForceRHR`? This answer actually seems alright. May 21, 2015 at 0:21
• @JakubKania I tried to come up with a `ST_ForceRHR` solution, but was not successful. Would be thankful for hints... I tried `ST_Dump(ST_Boundary(ST_ForceRHR(the_geom)))` Oct 17, 2015 at 11:00