Calculate the angle of exterior rings PostGIS (polygons & multipolygons)

postgis spatial functions
postgis polygon search
postgis box
postgis polygon to multipolygon
postgis point intersect polygon
postgis geometry vs polygon
st_union
postgis geography

SAMPLE DATA:

CREATE TABLE poly_and_multipoly (
  "id" SERIAL NOT NULL PRIMARY KEY,
  "name" char(1) NOT NULL,
  "the_geom" geometry NOT NULL
);
-- add data, A is a polygon, B is a multipolygon
INSERT INTO poly_and_multipoly (name, the_geom) VALUES (
    'A', 'POLYGON((7.7 3.8,7.7 5.8,9.0 5.8,7.7 3.8))'::geometry
    ), (
    'B',
    'MULTIPOLYGON(((0 0,4 0,4 4,0 4,0 0),(1 1,2 1,2 2,1 2,1 1)), ((-1 -1,-1 -2,-2 -2,-2 -1,-1 -1)))'::geometry
);

I have a table of multipolygons and polygons and I am trying to calculate the interior angles of the exterior rings in the table (i.e. no interior rings...) using ST_Azimuth. Is there any way to modify the attached query to use ST_Azimuth on the sp and ep of the linestrings?

    SELECT id, name, ST_AsText( ST_MakeLine(sp,ep) )
FROM
   -- extract the endpoints for every 2-point line segment for each linestring
   (SELECT id, name,
      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 id, name, (ST_Dump(ST_Boundary(the_geom))).geom
       FROM poly_and_multipoly
       ) AS linestrings
    ) AS segments;

1;"A";"LINESTRING(7.7 3.8,7.7 5.8)"
1;"A";"LINESTRING(7.7 5.8,9 5.8)"
1;"A";"LINESTRING(9 5.8,7.7 3.8)"
2;"B";"LINESTRING(0 0,4 0)"

The use of ST_Azimuth function at aengus example is giving an error because you can't use two generate_series as parameters for a single function. If you use ep / sp values in another subquery then you won't have any trouble.

Here is my solution:

-- 3.- Create segments from points and calculate azimuth for each line.
--     two calls of generate_series for a single function wont work (azimuth).
select id,
       name,
       polygon_num,
       point_order as vertex,
       --
       case when point_order = 1
         then last_value(ST_Astext(ST_Makeline(sp,ep))) over (partition by id, polygon_num)
         else lag(ST_Astext(ST_Makeline(sp,ep)),1) over (partition by id, polygon_num order by point_order)
       end ||' - '||ST_Astext(ST_Makeline(sp,ep)) as lines,
       --
       abs(abs(
       case when point_order = 1
         then last_value(degrees(ST_Azimuth(sp,ep))) over (partition by id, polygon_num)
         else lag(degrees(ST_Azimuth(sp,ep)),1) over (partition by id, polygon_num order by point_order)
       end - degrees(ST_Azimuth(sp,ep))) -180 ) as ang
from (-- 2.- extract the endpoints for every 2-point line segment for each linestring
      --     Group polygons from multipolygon
      select id,
             name,
             coalesce(path[1],0) as polygon_num,
             generate_series(1, ST_Npoints(geom)-1) as point_order,
             ST_Pointn(geom, generate_series(1, ST_Npoints(geom)-1)) as sp,
             ST_Pointn(geom, generate_series(2, ST_Npoints(geom)  )) as ep
      from ( -- 1.- Extract the individual linestrings and the Polygon number for later identification
             select id,
                    name,
                    (ST_Dump(ST_Boundary(the_geom))).geom as geom,
                    (ST_Dump(ST_Boundary(the_geom))).path as path -- To identify the polygon
              from poly_and_multipoly ) as pointlist ) as segments;

I've added some complexity as I wanted to identify each polygon to avoid mixing lines for angle calculation.

As sqlfiddle don't have support for PostGIS I've uploaded this example with some complementary code to github here

Chapter 14. PostGIS Special Functions Index, ST_Perimeter - Returns the length of the boundary of a polygonal geometry or ST_Azimuth - Returns the north-based azimuth as the angle in radians ST_Roughness - Returns a raster with the calculated "roughness" of a DEM. ST_ExteriorRing - Returns a LineString representing the exterior ring of a Polygon. geometry ST_ExteriorRing(geometry a_polygon); Description. Returns a line string representing the exterior ring of the POLYGON geometry. Return NULL if the geometry

You can add the azimuth calculation into your subquery as follows. Note that ST_Azimuth computes the angle clockwise from down to up, so more work will be required to ensure that az is actually the interior angle.

     SELECT id, name, ST_AsText( ST_MakeLine(sp,ep) ), az
    FROM
      -- extract the endpoints for every 2-point line segment for each      linestring
 (SELECT id, name,
  ST_PointN(geom, generate_series(1, ST_NPoints(geom)-1)) as sp,
  ST_PointN(geom, generate_series(2, ST_NPoints(geom)  )) as ep,
    ST_Azimuth(ST_PointN(geom, generate_series(1, ST_NPoints(geom)-1)),
    ST_PointN(geom, generate_series(2, ST_NPoints(geom)  ))) as az
FROM
   -- extract the individual linestrings
  (SELECT id, name, (ST_Dump(ST_Boundary(the_geom))).geom
   FROM poly_and_multipoly
   ) AS linestrings
) AS segments;

ST_ExteriorRing, Returns a line string representing the exterior ring of the POLYGON geometry. --If you have a table of polygons SELECT gid, ST_ExteriorRing(the_geom) AS  ST_ExteriorRing — Returns a line string representing the exterior ring of the POLYGON geometry. Return NULL if the geometry is not a polygon. Return NULL if the geometry is not a polygon. Will not work with MULTIPOLYGON

I needed the internal angles for only polygons. The accepted solution by vamaq did not consistently give me the internal angles, but sometimes external angles, for an arbitrarily shaped polygon. This is because it calculates the angle between two lines without any consideration of which side of those two lines is internal, so if the internal angle is obtuse vamaq's solution gives the external angle instead.

In my solution I use the cosine rule to determine the angle between two lines without reference to the Azimuth and then use ST_Contains to check whether the angle of the polygon is acute or obtuse so that I can ensure the angle is internal.

CREATE OR REPLACE FUNCTION is_internal(polygon geometry, p2 geometry, p3 geometry)
RETURNS boolean as
$$
    BEGIN
        return st_contains(polygon, st_makeline(p2, p3));
    END
$$ language plpgsql;


CREATE OR REPLACE FUNCTION angle(p1 geometry, p2 geometry, p3 geometry)
RETURNS float AS
$$
    DECLARE
        p12 float;
        p23 float;
        p13 float;
    BEGIN
        select st_distance(p1, p2) into p12;
        select st_distance(p1, p3) into p13;
        select st_distance(p2, p3) into p23;
        return acos((p12^2 + p13^2 - p23^2) / (2*p12*p13));
    END
$$ language plpgsql;

CREATE OR REPLACE FUNCTION internal_angle(polygon geometry, p1 geometry, p2 geometry, p3 geometry)
RETURNS float as
$$
    DECLARE
        ang float;
        is_intern boolean;
    BEGIN
        select angle(p1, p2, p3) into ang;
        select is_internal(polygon, p2, p3) into is_intern;
        IF not is_intern THEN
            select 6.28318530718 - ang into ang;
        END IF;
        return ang;
    END
$$ language plpgsql;


CREATE OR REPLACE FUNCTION corner_triplets(geom geometry)
RETURNS table(corner_number integer, p1 geometry, p2 geometry, p3 geometry) AS
$$
    DECLARE
        max_corner_number integer;
    BEGIN
        create temp table corners on commit drop as select path[2] as corner_number, t1.geom as point from (select (st_dumppoints($1)).*) as t1 where path[1] = 1;
        select max(corners.corner_number) into max_corner_number from corners;
        insert into corners (corner_number, point) select 0, point from corners where corners.corner_number = max_corner_number - 1;
        create temp table triplets on commit drop as select t1.corner_number, t1.point as p1, t2.point as p2,  t3.point as p3 from corners as t1, corners as t2, corners as t3 where t1.corner_number = t2.corner_number + 1 and t1.corner_number = t3.corner_number - 1;
        return QUERY TABLE triplets;
    END;
$$
LANGUAGE plpgsql;

CREATE OR REPLACE FUNCTION internal_angles(geom geometry)
RETURNS table(corner geometry, angle float)
AS $$
    BEGIN
        create temp table internal_angs on commit drop as select p1, internal_angle(geom, p1, p2, p3) from (select (c).* from (select corner_triplets(geom) as c) as t1) as t2;
        return QUERY TABLE internal_angs;
    END;
$$
LANGUAGE plpgsql;

Usage:

select (c).* into temp from (select internal_angles(geom) as c from my_table) as t;

Chapter 8. PostGIS Reference, ST_IsPolygonCCW — Tests if Polygons have exterior rings oriented counter-​clockwise These functions compute measurements of distance, area and angles. The angles of a regular polygon are equivalent, and their sides are as well. The sum of the exterior angles of a regular polygon will always equal 360 degrees. To find the value of a given exterior angle of a regular polygon, simply divide 360 by the number of sides or angles that the polygon has.

Chapter 7. PostGIS Reference, ST_InteriorRingN — Return the Nth interior linestring ring of the polygon geometry. ST_Azimuth — Returns the angle in radians from the horizontal of the vector defined by ST_Covers — Returns 1 (TRUE) if no point in Geometry B is outside Geometry Calculations are in the Spatial Reference System of this Geometry. The measure of each interior angle of a regular n-gon is. 1/n ⋅ (n - 2) ⋅ 180 ° or [(n - 2) ⋅ 180°] / n. The sum of the measures of the exterior angles of a convex polygon, one angle at each vertex is. 360 ° The measure of each exterior angle of a regular n-gon is. 360 ° / n

Chapter 6. PostGIS Reference, Calculations are in the Spatial Reference System of this Geometry. The optional third Return the exterior ring of the polygon geometry. Return Rotate the geometry around the Z, X or Y axis by the given angle given in radians. Follows the  The dominant angle of a polygon is the angle of longest collection of segments that have similar orientation. This angle will be stored in the specified field in decimal degrees from true north. Use this tool to determine the trend of a polygon and use the resulting angle to orient symbology such as markers or hatch lines within the polygon.

PostGIS 2.0 Cheat Sheet v2.0.2, PostgreSQL PostGIS Geometry/Geography/Box Types. Abstract ST_ExteriorRing — Returns a line string representing the exterior ring of the POLYGON geometry. ST_InteriorRingN — Return the Nth interior linestring ring of the polygon geometry. Calculations are in the Spatial Reference System of this Geometry. If you just want majority orientation, check out @Mapperz 's answer above. Otherwise, as you say you could split the polygons into lines using the Polygon To Line tool. This adds a left and right FID field, where the field is -1 if there is no outer polygon - this can cause a bit of mucking about though if your buildings are adjacent or overlap.

Comments
  • wow! This is very helpful (and very well understandable code logic). Thank you.
  • Thanks. hm. I'm getting an error with that...guess I'll have to debug... 'ERROR: functions and operators can take at most one set argument'