9

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.

enter image description here

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?
  • 2
    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) – Tamas Kosa May 18 '15 at 15:34
  • 1
    @TamasKosa: You have the essence of a good answer. Why not expand it into one? Also, for normals, azimuths need +/- pi/2. – Martin F May 18 '15 at 17:33
  • 1
    @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 ;) – n1000 May 18 '15 at 19:03
  • Just curious, why do you want normals ... sun exposure? – Martin F May 19 '15 at 0:46
  • 1
    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. – John Powell May 19 '15 at 6:32
9

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: enter image description here

  • 1
    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. – n1000 May 19 '15 at 10:19
  • 1
    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. – Tamas Kosa May 19 '15 at 11:01
  • 2
    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. – n1000 May 19 '15 at 11:46
  • 1
    What about ST_ForceRHR? This answer actually seems alright. – Jakub Kania May 21 '15 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))) – n1000 Oct 17 '15 at 11:00

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.